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ABSTRACT 


This report deals with progress made on the Grant NSG-3048 during the 
calendar year beginning M a r ch 1, 1977 and ending February 28, 1978. This 
year coincides with Supplement No. 2 of the award, which originated on 
March 1, 1975. The NASA Technical Officer for this period was Dr. Bruce 
Lehtinen of Lewis Research Center. The directors of the research at the 
University of Notre Dame were Dr. R. Jeffrey Leake and Dr. Michael K. Sain 

General goals of the research have been classified into two cate- 
gories. The first category involves the use of modem multivariable fre- 
quency domain methods for control of engine models in the neighborhood of 
a quiescent point. The second category involves the use of nonlinear mod- 
elling and optimization techniques for control of engine models over a 
more extensive part of the flight envelope. 

Substantial progress has been made in both categories. 

In the frequency domain category, works have been published in the 
areas of low- interaction design, polynomial design, the CARDIAD* method, 
and multiple setpoint studies. A number of these ideas have progressed 
to the point at which they are starting to attract practical interest. 
Further effort is yet required, however, to carry the ideas to maturity 

*The acronym stands for .Complex Acceptability Region for DlAg onal Dom- 
inance. See report for details. 


and to ensure their adequate dissemination. A highlight of the year was 
the incorporation of realistic jet engine data as a theme problem into 
the International Forum on Alternatives for Linear Multivariable Control. 
Dr. Sain was Program Chairman for this meeting, which attracted nearly 
two hundred persons from industry, laboratories, and universities to 
hear thirty papers focused in the general subject area of this grant. 

In the nonlinear category, advances have been made both in engine 
modelling and in the details associated with software for determination 
of time optimal controls. Nonlinear models for a two spool turbofan 
engine have been expanded and refined; and a promising new approach to 
automatic model generation has teen placed under study. A two time scale 
scheme has been developed to do two-dimensional dynamic programming, and 
an outward spiral sweep technique has greatly speeded convergence times 
in time optimal calculations. 

The details of these and other aspects of the year T s investigations 
may be found in the body of the report, which covers the most active 
grant period to date. 
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I. INTRODUCTION 


The purpose of this section is to provide some of the broad back- 
ground which underlies and clarifies the general nature or the Research 
Highlights, which are stated in the section following. 

Initiation of Grant NSG-3048 in March 1975 was timed with develop- 
ments in the engine industry, which was beginning to experience some lim- 
itations in the application of classical hydromechanical control tech- 
nique as the primary base technology for modern engines with ever in- 
creasing sophistication. At the same time, milestone developments in 
digital hardware began to open realistic possibilities for onboard comp- 
utation to an extent not heretofore possible. This confluence of events 
led directly to the concept of increasing the role of electronics in 
engine control. In turn, the availability of digital electronics itself 
created a wide variety of opportunity for application of new control de- 
sign philosophy and technique. Among the earliest of such studies is 
the F100 Multivariable Control Synthesis Program [1] sponsored by the 
National Aeronautics and Space Administration, Lewis Research Center 
and the Air Force Aero-Propulsion Laboratory, Wright-Patterson Air Force 
Base. This program is currently in the test phase. 

The advent of digital technology on the engine scene offers not 
only the opportunity to control more engine Variables but also the pos- 
sibility of integrating engine and airframe control. Studies of this 
type have also begun. 

Primary tools in the F100 Multivariable Control Synthesis Program 


were linear quadratic regulator (LQR) theory in the linear case.. For the 
global control, nonlinear optimal methods were not directly applied. 

The purpose of Grant NSG— 3048 is to evaluate alternatives to LQR in 
the linear case and to examine nonlinear modelling and optimization ap- 
proaches for global control. 

Context for the studies is set by the DYNGEN digital simulator [2] . 
Based upon earlier computer codes GENENG [3] and GENENG II [4], DYNGEN 
has the combined capabilities of [3] and [4], for calculating steady- 
state performance, together with the further capability for calculating 
transient performance. DYNGEN uses a modified Euler method to solve the 
differential equations which model the dynamics of the engine. This mod- 
ified Euler method permits the user to specify large time steps, for ex- 
ample a tenth of a second; and this can result in considerable savings of 
execution time. On the other hand, convergence problems are sometimes 
encountered with DYNGEN wher small time steps are used. 

The DYNGEN digital simulation is particularized to a given situation 

by a process of loading data for the various maps associated with a given 

engine. The maps for the Grant NSG-3048 have been provided by engineer- 

» 

ing personnel at Lewis Research Center. These maps correspond to- a 
paper engine, which is not closely identified with any current engine. 

But the data do correspond in a broad, general sense to realistic two 
spool turbofan engines. The simulation provides for two essential con- 
trols, main burner fuel flow and jet exhaust area. Portions of the en- 
velope which can be used for linear or nonlinear experimentation are 
limited by the convergence capabilities of the available engine data on 
DYNGEN. 
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With respect to multivariable frequency domain work, the basic ap- 
proaches may be classified into two groups. These two groups are often 
called "direct" and "indirect". 

* The direct approach can usually be recognized by its attention to 

achieving a completely specified dynamic performance. Such ideas have 
been discussed from the early days of organized control study. See, 
for example, [5] and [6]. In fact, some of the earliest attempts to ex- 
pand the direct approach to the multi- input, multi-output case involved 
work with jet engines [7,8 J. Direct approaches in multivariable appli- 
cations typically involve matrices of transfer functions. In the 1950 T s, 
there were some nontrivial difficulties with such methods in cases of 
more than one input and output. Among these difficulties may he mentioned 

(1) the meaning and extent of cancellations of various types 
in the transfer functions, 

(2) the question of loop stability, 

(3) the problem of specification, and 

(4) numerical computations. 

Advances in the last two decades have resolved (1) and (2); in industry, 
a reservoir of expertise has built up relative to (3); and progress in 
the numerical sciences is rapidly achieving an interface witb. control 
theory in such a sany as to resolve (4) . It is believed that work on 
this grant is continuing to stimulate rapid developments in areas ; (3^ and 

Indirect approaches are usually recognizable by tha$f affinities 


to the classic worlcs of Nyquist or Evans, involving, respectively, fre- 
quency response plats in the complex plane or versions of the root locus. 
Xt is a relatively easy matter to describe the focus of generalized 
, Nyquist methods. The key constituent ideas are related to three poly- 
nomials : 

(1) . p (s) - the closed loop characteristic polynomial (CLCP) , 

c 

(2) p (s) - the open loop characteristic polynomial (OLCP) , and 

o 

(3) jM(s) [— the determinant of return difference. 

The CLCP is a polynomial whose zeros characterize the exponential im- 
pulse response of the closed loop control system; the OLCP serves the same 
purpose for the open loop system. M(s) is a matrix of transfer functions 
associated with the following experiment. Break the control loops at a 
convenient point and inject impulses. The difference between the trans- 
form of the signal injected and that which returns at the other end of 
the loop is established by the columns of M(s) . The quantities p^Cs), 
P o (s) and M(s) derive their importance from the fact that they are re- 
lated to each other by the equation 

P c (s) = jM(s) | P D (s). 

Typically, P c (s) is known; and M(s) is partly given and partly designed, 
in such a way that P c (s) becomes desirable. 

Generally speaking, a Nyquist plot of |m(s) ] tends to contain the 
same types oT information which proved so successful in classical designs. 
A great deal of the design effort centers upon the way in which dynamical 
compensation affects the determinant which acts on M(s) . There are three 
well recognized ways to study this effect. These are 
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(1) direct construction of |M(s) j by any of the known methods 
for determinant calculation, 

(2) construction of the eigenvalues of Jm(s) j as a function of 
s, and use of the idea that the determinant is equal to 
the product of its eigenvalues [9], and 

(3) ‘ design of compensation so that M(s) is approximately dia- 

gonal, and establishment of a relation between the plot 
of |M(s) [ and plots of the diagonal elements of M(s) [10]. 

It is believed that work on this grant has advanced the application of 
all these methods to jet engine design, but particularly method (3), 
where a special technique has been developed to design compensation so 
that M(s) is approximately diagonal. This technique is called the CAEDIAD 
Plot, where the acronym stands for Complex Acceptability Region for D1A- 
gonal Dominance , the latter term referring to a specific definition of 
"approximately diagonal." 

With respect to nonlinear modelling and optimization, the emphasis 
has been twofold: to develop good analytical nonlinear models of the 
jet engine and to use these models in conjunction with techniques of 
mathematical programming in order to develop advances in global control 
over significant reaches of the flight envelope. 

In general, there are several aspects to this part of the investi- 
gation. First, it is possible to conceive the basic differential equa- 
tions from fundamental principles. In this case, there are usually about 
sixteen nonlinear differential equations, as well as a large number of 
nonlinear static functions which serve as part of the coupling between 


the equations. These functions often have more than one argument. If 
the equations arise in this fashion, then there is a significant need 
to identify the parameters. This must normally be done from the DYNGEN 
digital simulation. Second, it is possible to assume a general form for 
the nonlinear differential equations in such a way that fundamental prin- 
ciples are not ignored but that added emphasis is placed upon general 
mathematical form. If this general form is chosen according to a scheme 
designed to make Maximum use of the type of data which is directly avail- 
able from the digital simulation, then a type of "automatic" nonlinear 
model generation becomes possible. Third, whether the first or second 
modelling procedure is employed, there is almost always a need to con- 
sider the problem of reducing the order of the models. Though order 
reduction can often he highly mathematical in nature, it is almost al- 
ways the case that the reduced order model depends upon the scaling Of 
the equations. As a result, the final reduced models often depend in a 
nontrivial way upon physical insight, as well as mathematical method. 

Work on this grant has focused especially upon the first and second 
aspects of the modelling problem, with a gradual specialization toward 
automatic model generation. 

Insofar as optimization is concerned, the stress has been placed 
upon time optimal control, and considerable effort has been invested in 
specialized programming methodology designed to take maximum advantage 
of the particular features of jet engine models. , 

In the next section, the highlights of activities carried out during 
the calendar year corresponding to this report are presented. 
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II. HIGHLIGHTS OF THE RESEARCH 


This section is a brief statement of the main achievements under 
, Grant NSG-3048 during the period from March 1, 1977 to February 28, 1978. 
There are too major subdivisions, according to the main thrusts of the 
investigation. The first of these is Local Multivariable Frequency Domain 
Methods; and the second is Global Nonlinear Optimal Methods. 

For the most part, the wording of these paragraphs has been con- 
strained so as to be as nontechnical as possible. Nonetheless, some 
readers may find it useful to review the basic introduction provided in 
Section I. 

A. Local Multivariable Frequency Domain Methods 

During the calendar year ending on February 28, the following results 
were achieved in the area of modern, frequency domain control of turbofan 
engine models. 

Cl) The first formal documentation of the CARDIAD method 
(Complex Acceptability Region for DIAg onal Dominance) 
was completed. See (1), Section III. Though supported 
principally under a theory grant from the National 
Science Foundation, this technique had its origin in 
class studies of older methods for approximate de- 
coupling of jet engine models in the frequency domain. 

Almost all of the examples in this thesis were taken 
from FlOO-like engine data. 

(2) The first documented studies of direct Nyquist plots 
of return difference determinants for jet engine 
models were completed. See (3), Section III. This 
thesis has been a helpful ancillary tool in general 
frequency domain design. 

(3) The first frequency domain closed loop compensation 
and simulation of a DYNGEN turbofan engine model was 
achieved. See (5), Section III. As explained in 
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Section I, the DYNGEN simulation supplies two control 
inputs . 

(4) The first study of polynomial techniques for exact 
model matching control of jet engine models in the 
frequency domain was reported. See (6), Section III. 
This paper has been pivotal in promoting the numer- 
ical advance of such techniques for applications. 

More will be reported in the subsequent semi-annual 
status report. 

(5) About a given design point, linear models of the 
standard type are obtained from the DYNGEN simulation 
by the DYGABCD routine [11] . In order to use DYGABCD 

at off-design points, however, modifications to DYGABCD 
necessary to the research had to be accomplished. See 
(7) , Section III. 

(6) The CARDIAD plot was applied to a series of DYNGEN off- 
design point models in order to determine its utility 
as a method for global classification of interaction 
characteristics of jet engines. See (14), Section III. 
The results were positive. 

(7) An entire conference was convened from industry, lab- 
oratories, and universities to hear speakers from sev- 
eral countries apply their theories to a theme problem 
developed from jet engine data. See (15), Section III. 
This meeting resulted in a book publication [12] . 

(8) The CARDIAD methodology was extended to the three-con- 
troi-input case and applied successfully to Pratt- 
Whitney data for the F100 engine. See (16) , Section 
III. 

(9) A joint seminar series was established between the 
Department of Electrical Engineering at Notre Dame and 
the Energy Controls Division of the Bendix Corporation 
at South Bend, in areas of mutual interest. This has 
resulted in published work. See (17), Section III. 


B. Global Nonlinear Optimal Methods 

The major advances and results achieved during the past year in the 


area of global nonlinear optimal methods are the following. 


(1) The hierarchy of analytical nonlinear models for the 
two spool turbofan jet engine has been expanded and 
refined. See (2), Section III. This effort is in 
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keeping x*ith general interest in the industry concern- 
ing improvement of compact general models. 

(2) A comparison has been achieved between the use of a 
linear affine model and a nonlinear model for time 
optimal control studies of a single spool engine. See 
(4), Section III. The results again support the search 
for reasonably simple nonlinear models, in the sense 
that they argue in favor of models whose nonlinearity 
is not excessively complicated. 

(3) ' A method using linear quadratic regulator methods to 

obtain decoupled control has been tested on various 
engine models. See (8), Section III. This is also 
an outgrowth of the joint Notre Dame - Bendix seminar 
series mentioned in (9), Section IIA, above. 

(4) A two time-scale scheme has been developed in order to 
do two-dimensional dynamic programming on a fifth 
order model of a jet engine. See (10), Section III. 

This is part of a continuing study of time-optimal 
control methods applied to nonlinear engine models. 

(5) Convergence times in time-optimal successive approx- 
imation dynamic programming have been dramatically 
improved through development of a scheme for a spiral 
out sweep from the target. See (10), Section III. 

(6) A completely automatic method for obtaining nonlinear 
analytical models for engine simulations has been de- 
veloped and tested numerically. See (11), Section III. 

This approach offers considerable promise for improve- 
ment over previous methods. 

(7) A discrete maximum principle has been developed for 
nonlinear systems having the property that the con- 
trol is constrained by the present state. See (12), 

Section III. 

(8) A family of optimal feedback control laws has been 
developed and simulated for a variety of models. 

Further details concerning these highlights may be found in Sections IV 
and V. Also, as described in Section III following, a number of the doc- 
uments have been included as appendices. 

The next section contains a list of publications completed during 
the current year of the grant. 
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III. PUBLICATIONS 


This section provides a list of the nineteen documents completed 
» during the year March 1, 1977 through February 28, 1978. The works are 
ordered chronologically. 

Some of the listings are followed by an alphabetical code consisting 
of one or more of the letters A, M, and R. The letter A signifies that 
the document or an abstract thereof appears as one of the appendices to 
this report; the letter M signifies that the document comprises a thesis 
for the degree of Master of Science in Electrical Engineering; and the 
letter R declares that the document summarizes an effort which was closely 
integrated with, but not directly supported by, the activities of this 
grant. 


Completed publications from earlier years are not included in this 
list; but a total listing of all the grant documents has been provided 
as an appendix to the report. See the Table of Contents. 

(1) R.M. Schafer, "A Graphical Approach to System Dominance”, 
Technical Report EE 772, University of Notre Dame, April 
1, 1977. (M,R) 

(2) W.E. Longenbaker and R.J. Leake, ’'Hierarchy of Simula- 

tion Models for a Turbofan Gas Engine”, Proceedings 
Eighth Annual Pittsburgh Conference on Modeling and 
Simulation, April 1977. (A) 

(3) P.W. Hoppner, "The Direct Approach to Compensation of 
Multivariable Jet Engine Models”, Technical Report EE 
774, University of Notre Dame, May 1977, (M,R) 

(4) R.R. Gejji and R.J. Leake, "Time-Optimal Control of a 
Single Spool Turbojet Engine Using a Linear Affine 
Model”, Technical Report EE 7711, University of Notre 
Dame, June 1977. 


10 


(5) R.M. Schafer, R.R. Gejji, P.W. Hoppner, W.E. Longenbakei 
and M.K. Sain, "Frequency Domain Compensation of a DYNGEN 
Turbofan Engine Model", Proceedings Sixteenth Joint Auto- 
matic Control Conference , pp. 1013-1018, June 1977. (A) 

(6) R.R. Gejji and M.K. Sain, "Application of Polynomial 
Techniques to Multivariable Control of Jet Engines", Pro- 
ceedings Fourth IFAC Symposium on Multivariable Tech- 
nological Systems , pp. 421-429, July 1977. (A) 

(7) R.R. Gejji, "Use of DYGABCD Program at Off-Design Points", 
Technical Report EE 7703, University of Notre Dame, July 
1977. 

(8) E.A. Sheridan and R.J. Leake, "Non- Interactive State Re- 
quest Jet Engine Control with Non-Singular B Matrix", 
Proceedings Twentieth Midwest Symposium on Circuits and 
Systems , pp. 539-543, August 1977. (R) 

(9) R. Gejji, R.M. Schafer, M.K. Sain, and P. Hoppner, "A Com- 
parison of Frequency Domain Techniques for Jet Engine Con- 
trol System Design", Proceedings Twentieth Midwest Sympo- 
sium on Circuits and Systems , pp. 680-685, August 1977. (A) 

(10) W.E. Longenbaker and R.J, Leake, "Time Optimal Control of 

a Tero-Spool Turbofau Jet Engine", Technical Report EE 7714, 
University of Notre Dame, September 1977. (A,M) 

(11) R.J. Leake and J.G. Comiskey, "A Direct Method for Obtain- 
ing Nonlinear Analytical Models of a Jet Engine", Proceed- 
ings International Forum on Alternatives for Linear Multi- 
variable Control , National Electronics Conference, Chicago, 
pp. 203-212, October 1977. (A) 

(12) J.A. Ortega and R.J. Leake, "Discrete Maximum Principle 
with State Constrained Control”, SIAM Journal on Control 
and Optimization , Vol. 15, No. 6, pp. 109-115, November 
1977. (R) 

(13) Michael K. Sain and V. Seshadri, "Pole Assignment and a 
Theorem from Exterior Algebra", Proceedings IEEE Conference 
on Decision and Control , pp. 291-295, December 1977. (R) 

(14) R. Michael Schafer and Michael K. Sain, "Some Features of 
CARDIAD Plots for System Dominance", Proceedings IEEE Con- 
ference on Decision and Control , pp. 801-806, December 

1977. (A) 

(15) M.K. Sain, "The Theme Problem", in Alternatives for Linear 
Multivariable Control , M.K. Sain, J.L. Peczkowski and J.L, 
Melsa, Editors. Chicago: National Engineering Consortium, 

1978, pp. 20-30. (A) 
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(16) R.M. Schafer and M.K. Sain, "Input Compensation for 
Dominance of Turbofan Models", in Alternatives for 
Linear Multivariable Control , M.K. Sain, J.L. 

Peczkowski, and J.L. Melsa, Editors. Chicago: Na- 
tional Engineering Consortium, 1978, pp. 156-169. (A) 

(17) J.L. Peczkowski and M.K. Sain, "Linear Multivariable 
Synthesis with Transfer Functions", in Alternatives 
for Linear Multivariable Control , M.K. Sain, J.L. 
Peczkowski, and J.L. Melsa, Editors. Chicago: Na- 
tional Engineering Consortium, 1978, pp. 71-87. (R) 

(18) R.J. Leake and M.K. Sain, "Semi-Annual Status Report, 
NASA Grant NSG-3048, 'Alternatives for Jet Engine Con- 
trol', Supplement No. 2", March 1, 1977-August 31, 1977. 

(19) R.J. Leake and M.K. Sain, "Final Technical Report, NASA 
Grant NSG-3048, 'Alternatives for Jet Engine Control', 
Supplement No. 2", March 1, 1977-February 28, 1978. 


IV. LOCAL MULTIVARIABLE FREQUENCY DOMAIN METHODS 

Progress on local multivariable frequency domain methods has been 
achieved during this grant period in the areas of Low Interaction Design, 
Polynomial Design, Extension of CARDXAD Method, and Multiple Setpoint 
studies. 

Low Interaction Design 

As mentioned in the Final Report for NASA Grant NSG-3048, Supplement 
No. 1, a promising new technique for designing dynamical compensation be- 
gan to develop in the Fall of 1976. This methodology, built upon what 
are currently being called CARDIAD plots, was only being tentatively con- 
sidered in October, 1976 when the continuation proposal for NASA Grant 
NSG-3048, Supplement No. 2, was being written. Based upon favorable pre- 
liminary reaction by personnel from NASA Lewis Research Center, a decision 
was made to investigate further the use of CARDIAD plots as a design aid 
for turbofan engine control in the frequency domain. In essence, this 
study proved to be successful enough that it really dominated the re- 
maining time period of Supplement No. 1 and has continued through Sup- 
plement No. 2. 

A great deal of the power of the CARDIAD plot arises from its sim- 
plicity. For each frequency, a circle is constructed on a planar plot. 
Data for the center and radius of this circle is obtained from the com- 
plex transfer function matrix of the plant. The circle may be solid or 
dashed. If solid, the inside of the circle defines the acceptable com- 
plex region for the value of a frequency dependent compensator element 
in order to achieve dominance. If dashed, the outside of the circle de- 
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fines the acceptable region. As the frequency follows a standard Nyquist 
pattern, these circles result in a CARDIAD plot. (Complex Acceptability 
Region for DIA gonal Dominance)* This plot has been shown to speak con- 
structively to the issue of compensator choice to reduce interaction. 

As an example of the CARDIAD plot application to the turbofan engine 
control .problem, a linear model obtained using DYGABCD on the DYNGEN dig- 
ital engine simulator was used to illustrate control design at the 1977 
Joint Automatic Control Conference. The paper based upon this effort, 
which may be seen in Appendix C, utilizes a two-input, five-state, two- 
output engine model in which the inputs are fuel flow and nozzle area, 
the states are compressor speed, fan speed, burner exit pressure, after- 
burner exit pressure, and high turbine inlet energy, and the outputs are 

.thrust and high turbine inlet temperature. 

r 

Typical examples of CARDIAD plots for such engine models may be seen 
in Figures 2-5 of Appendix C. The investigators involved in this study 
have seen the same type of plots arising from a variety of engine data. 
This has raised the interesting question of whether there may be a mean- 
ingful concept of "engine interaction footprint" in the sense of the 
CARDIAD plot. 

Of particular interest is the plot shown in Figure 3 of Appendix C. 
Students of classical control theory will immediately recognize the near 
semicircular nature of this plot. Such semicircular behavior has been 
observed frequently and serves to specify a sort of essential lead-lag 
classical compensator* element which can achieve diagonal dominance. 

Using the CARDIAD approach, it has been possible to achieve diagonal 
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dominance at all frequencies on typical engine models. Moreover, only 
simple compensators have been required to do this. While it is not pos- 
sible to apply the same degree of credibility to the model itself at all 
frequencies, it is nonetheless of considerable theoretical interest to 
be able to make this accomplishment, especially for (A, B, C, D) type 
plant models which have the D matrix present to approximate modelling 
errors at high frequencies. Further insight into the significance of 
these steps can be obtained by examining Appendix B of the Final Report 
for Supplement No. 1. 

Appendix C of this report also contains evidence of two other 
facets of the applications researches conducted under this grant. Fig- 
ures 6-8 are characteristic locus plots for the plant, after it was com- 
pensated by the CARDIAD methodology. This combination of ideas, namely, 
the CARDIAD plot and the characteristic locus, has been quite helpful! 
in studies conducted up to this time. Softwares were developed and ex- 
perience gained with the characteristic locus on the original grant NSG- 
3048, as well as on Supplement No, 1. 

Of additional interest also are Figures 9-11 of Appendix C. These 
figures deal with an aspect of frequency domain control research which 
may he called a "direct” approach. The term "direct" refers to a direct 
construction of the determinant of return difference as frequencies fol- 
low a standard Nyquist pattern. This tool really underlies all the mod- 
ern frequency domain ideas; but it is usually handled obliquely, as for 
example by the CARDIAD idea or by the characteristic locus. Studies of 
the direct approach to determinant of return difference have confirmed 
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that it is at the very least a revealing analysis method. Figure 11 of 
Appendix C, for example, reveals a condition of gain limitation which is 
understandable in a global way not so easily visualized by separate char- 
acteristic locus plots. It should be emphasized, moreover, that diagrams 
of the type of Figure 11 can be drawn without any regard of plant size 
in terms of inputs, states, and outputs. 

Efforts to use the direct approach for design, as well as for anal- 
ysis, have posed nontrivial algebraic questions. Some insight has been 
gained, but no breakthroughs have occurred as yet. 

Polynomial Design 

The principal efforts and results obtained in applying polynomial 
design techniques to the turbofan engine control problem have been re- 
ported In the Final Report on Supplement No. 1. 

During the present grant period, a paper on this work was presented 
at the International Federation of Automatic Control's Fourth Symposium 
on tfultivariable Technological Systems at the University of New Brunswick 
in Fredericton. See Appendix D. 

There is now little doubt that the control area is experiencing a 
resurgence of interest in transfer function methods. As part of this 
resurgence, grant work on polynomial design has pointed out the necessity 
of increased attention to numerical method. The investigators also be- 
lieve that it has stimulated other workers to begin numerical studies. 


The transfer function has a number of key properties which have long 
made it popular with control practitioners. For example, the transfer 
function is unique relative to similarity transformations on the state 



space 


Much work remains to be done, however, on computational aspects of 
transfer function design. 

A presentation comparing the design experiences of the investigators 
under this grant, in the frequency domain, was made at the 1977 Midwest 
Symposium on Circuits and Systems, Lubbock, Texas, in August. A copy of 
this brief manuscript may be seen in Appendix E. 

Extension of CARDIAL Method 

All the work so far mentioned in regard to the CARDIAL plot was 
carried out for plants having two inputs and two outputs. In this sit- 
uation, it is certainly true that the plots have many interesting prop- 
erties . 

On the basis of this e> lerience, which was in its final stages in 
May, 1977, it was decided to extend the CARDIAD theory to the three-input, 
three-output case. Also, the National Engineering Consortium’s Inter- 
national For um on Alternatives for Linear Multivariable Control, which 
took place in October, 1977, offered a prime opportunity to apply the 
theoretical extension, inasmuch as the Forum contains a Theme. Problem 
based on a linearized model of a modern tur'bofan engine. See Appendix 
I. Assistance in the theoretical extension of CARDIAD to the three-in- 
put, three-. mtput case was provided by support extended to Dr. Sain by 
the National Science Foundation under Grant ENG 75-22322. NASA support 
under Supplement No. 2 was focused on the turbofan application. 

Appendix J, "Input Compensation for Dominance of Turbofan Models", 
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provides a complete description of the successful work on this extension 
and its application. 

Technically, the extension of CARD I AD plot methods to the three-in- 
put, three-output case involved the use of a bound which provides suf- 
ficient conditions for diagonal dominance. Many possibilities exist for 
the selection of such a bound, and there remains considerable opportunity 
for further research along these lines. After examination of a number 
of basic bounding possibilities, an initial selection was made in such a 
way that the bound will extend to the case of a plant having p inputs 
and p outputs, where p is any positive integer greater than or equal to 
two, and that the bound will be tight in the place, where it matters the 
most where the plant is close to failing the dominance test. 

From an engineering point of view, it was necessary to develop soft- 
ware to extend the CA5DIAD idea and to establish viewpoints for studying 
the plots in more complicated cases. The CARD I AD analysis was divided 
into two phases. The first phase assumed one off-diagonal compensator 
element to he zero. An advantage of such a phase lies in its conceptual 
reduction to the situation of Appendix C, where greater design experience 
is available. The first-phase approach was adequate for about half of 
the application to the National Engineering Consortium Theme Problem. 

The second phase assumed both off-diagonal compensator elements to be 
nonzero, but drew the plot in such a way that the designer could as- 
certain what would happen when one of those elements was zero. This sec- 
ond phase approach was successful in completing the design for the Theme 
Problem. 
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A number of noteworthy points should be brought out concerning this 

extension research. • . 

(1) These CARDIAD methods have been successful without assuming a fixed 
form of the compensator . The importance of this fact can scarcely 
be overemphasized. If a fixed form is assumed, it' may happen that 
the -form is inadequate to fit the essential plant characteristics; 
and, as a result, it can well be the case that essential insight 

is lost. 

(2) The CARDIAD approach, applicable to design localized to just one 

side of the plant, as for example the input, can be used to affect 
outputs that are not measurable. Other methods that use compensation 
both at plant input and plant output often depend upon moving the 
output compensator around the loop an operation which is not pos- 

sible unless those outputs drive the loop. It would seem that this 
could he quite important- in the case of key outputs such as high 
turbine inlet temperature and thrust. 

(3) In cases studied so far, and there have been over a dozen of them, 
the CARDIAD plot has achieved dominance over all frequencies, even 
when the plant has (A, B, C, D) form and transmission does not roll 
off to zero at high frequencies. This design power has been ac- 
companied by a need for only relatively simple dynamical compensation. 

(4) In practice, gain selection in compensators has to be done with some 
care, so as not to invalidate the accuracy of the linearized model. 

The CARDIAD approach conveys considerable direct insight into the 
gains available; and does not leave the choice indirectly to an op- 
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Because of these features, the investigators feel that the CARDIAD plot 
is helping to push hack the research frontier in frequency domain ap- 
proaches to approximate decoupling. 

Multiple Setpoints 

From the outset, the CARDIAD plot has offered much promise for the 
general control problem which involves linearization at multiple set- 
points, design at each setpoint, and a piecing together of these designs 
to achieve global effects. Such technique is certainly the norm both in 
present-day practice and in current research for the turbofan engine. 

Basically, the idea is to construct CARDIAD plots to each setpoint 
and to use these to study the interaction features of the nonlinear en- 
gine model over a more global operating regime. The investigators be- 
lieve that such studies can be helpful in selecting setpoints for design 
and in constructing compensation which works toward globa l dominance . 

Work has been proceeding along these lines, with the aid of set- 
points involving two inputs and two outputs from NASA's DYNGEN engine 
simulator. The first documentation can be found in Appendix 5. 

¥ 

In Appendix H, the setpoints are determined by fuel flows of 2.145, 
2.31, 2.475, 2,64, and 2.75 LBM/SEC. Figures 1-10 of the appendix con- 
tain the corresponding CARDIAD plots. Consider Figures 1, 3, 5, 7 and 9, 
which focus on the first column. Note their clear similarity. Next con- 
sider Figures 2, 4, 5, 8 and 10. Again note their clear similarity. Be- 
cause of this similarity, it was possible to design one simple compensa- 
tor to achieve diagonal dominance at all five setpoints. This compensa- 
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V. GLOBAL NONLINEAR OPTIMAL METHODS 


This section is concerned with some of the details of the third year of 
effort on the global nonlinear optimal part of the research. As in the 
previous year, this part is primarily concerned with the control of a 
taro-spool engine. 

There are three main aspects of the work: 

DYNGEN Simulator Operation 
Turbine Engine Modelling 
Nonlinear Optimal Control 

DYNGEN Simulator Operations 

The DYNGEN simulator, equipped with DYGABCD, has been useful in 
nearly all studies related to the grant, as it provides a "real world" 
testing ground for the various control methods under investigation. How- 
ever, it is costly and has limitations. Two such limitations are the fact 
that it is difficult to get convergence at low rotor speeds, and that only 
two controls (WEB and Ag) are readily available to the user. 

The automatic generation of AECD matrices enabled by DYGABCD has 
been invaluable. It is felt that our work has contributed to the overall 
development of the simulator through feedback provided, as for example, 
the simple modification suggested in Mr. Gejji’s memorandum. See (7), 
Section III . 

In the early stages of the work, DYNGEN was of primary concern, but 
it is now fairly routine and attention has turned to other areas. 


Turbine Engine Modelling 

This phase of the work began with analog computer studies of a single 
spool engine. Then a considerable effort was spent to obtain a good ana- 
lytical model for a two sp'-ol model using fundamental physical consider- 
ations. This study resulted in a hierarchy of models as reported in W.E. 
Longenbaker ’ s M.S. Thesis, Appendix F, and the paper by Longenbaker and 
Leake presented at the Pittsburgh Conference on Modelling and Simulation, 
Appendix B. As indicated in the results of Longeribaker’s Thesis, the mod- 
els obtained were disappointing. Even linear affine models appeared to 
fare better. See (4), Section XII. 

As a result of i Ills experience, the main emphasis in the work has 
now turned toward automatic generation of models by computer methods. The 
first effort in this direction is reported in the Chicago International 
Forum paper by Leake and Comiskey, Appendix G. The basic approach is to 
use an approximation of 

x - f(x,u) 

which is of the form 

/* 

x = A(x,u) (x - g(u)). 

This form seems to work very well for jet engine models. In the first 
place, there is always a unique equilibrium point for a given fixed con- 
trol u, so 

X = g(u) 


is the equilibrium equation which is all important for steady state anal- 
ysis. In the second place, jet engine A matrices rarely have poles 
(eigenvalues) at the origin and hence they are invertible. 
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This is a great help in computing such models. Suppose then that 


f(x,u) *-= A(x,u) (x-g(u)). 


and let 


(x , u ) 
e e 


be an equilibrium pair satisfying 


x e = g(u e ). 


Then 
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and 


3f 


„ (x , a ) « - A(x , u ) (u ) . 
3u e* a e* e 3u e 
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How it is well known that — corresponds to the approximate system A 

1 dX 


matrix in the steady state, so in our model, A(x,u) is a running approx- 
imation of the system A matrix and it can thus be approximated by meas- 
uring the A matrix at equilibrium points of interest and interpolating. 


3f 


A key point is, however, that — corresponds to the system B matrix 


in the steady state, and hence 


B( V V = " A( V V If ^e 5 


or 


3g 


.- 1 , 


a («_) = ~ A (x , u ) B(x , u ) . 
du e e e e e 


This is where the invertibility of A comes in. 


Thus, if there is an automatic method of finding A and B matrices 
(as we have in DYGABCD) then we have an easy way to get measurements of 
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A(x,u), g(u), and (u) . 

The term is very important because it is tbe DC gain of the linear 
model from control to state. To see this, consider the transfer func- 


tion relation 


X(s) = (sI-A) B U(s) 


Then the s = 0 DC gain relation is 

- A _1 B. 

It follows from the above discussion that use of the model form pre- 
scribed permits one to key in on 

Authentic Equilibrium Values 
Authentic A Matrix Values 
Authentic DC Gain Values 


for a global nonlinear model by measuring only 


Equilibrium Values 


A Ma trices 


B Matrices, 


These measurements can usually he made by automated methods. 

?!he Chicago paper, Appendix G, was a first attempt to use the ap- 
proach. J.G. Comiskey's M.S. Thesis is to use Hermite polynomials which 
are well matched to derivative requirements. 


Honlinear Optimal Control 

It is felt that W.E. Longenbaker did a comprehensive job of refin- 
ing our basic successive approximation Dynamic Programming scheme and 


25 


applying it to the models he studied. Details can he found in Tech- 
nical Report No. EE-7714 t Appendix F, and in our Semi-Annual Status 
Report for the period March 1, 1977 - August 31, 1977. 


VI. SPECIAL INITIATIVES RELATED TO GRANT WORK 


Two special initiatives were carried out during this year. The first 
was a special session at the 1977 Joint Automatic Control Conference, and 
the second was an entire meeting, the International Forum on Alternatives 
for Linear Multivariable Control. 

Joint Automatic Control Conference 

A session "Turbofan Engine Control" was put together for this con- 
ference. Co-Chairmen and Organizers were Drs. Michael K. Sain and H. 
Austin Spang. The papers are listed below. 


1. System Identification Principles Applied to Multivariable Control 
Synthesis of the F100 Turbofan Engine 

R.L. DeHoff and W.E. Hall, Jr. 

Systems Control, Inc. (Vt.) 

2. Failure Detection and Correction for Turbofan Engines 
H.A. Spang, III and R.C. Corley 

General Electric Company 

3. Frequency Domain Compensation of a DYNGEN Turbofan Engine Model 

R.M. Schafer, R.R. Gejji, P.W. Hoppner, W.E. Longenbaker, and M.K. Sain 
University of Notre Dame 

See Appendix C. 

4. The Application of the Routh Approximation Method to Turbofan Engine 
Models 

W. Merrill 

NASA Lewis Research Center 

5. Minimum-Time Acceleration of Aircraft Turbofan Engines 

F. Teren 

NASA Lewis Research Center 

6. Optimal Controls for an Advanced Turbofan Engine 

G. L. Slater 

University of Cincinnati 
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International Forum on Alternatives for Linear Multivariable Control 


In October, 1977, Dr. Sain was Program Chairman for an entire meeting 
focused in the general subject area of this grant. Approximately two hur- 
, dred persons attended from industry, laboratories, and universities. About 
thirty papers were presented, many by invited authorities of international 
stature. Nearly two-thirds of these addressed themselves to a Theme Prob- 
lem, Appendix I, which was derived from researches on this grant. Two 
publications resulted. The Proceedings contained contributed papers and 
abstracts of invited papers. The book [12] contained invited papers and 
those contributed papers which best fit in with the Forum Theme . 

The Forum Program appears on the next two pages. More information 
can be found in [12]. 
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Thursday, October 13 

SESSION 2, Origins of the Theme Problem 
Chairman: M. K. Sain, University of Notre Dame 


F 

£ 


{; 

f; 


**- 

l 


8:00 Engine Criteria and Models for Multivariable Control System Design 

&. D* Hackney and R. J* Miller, Pratt-Whitney Aircraft Group, and L* L. Small, Air Force Aero-Propulsion 
Laboratory 

8:30 A Practical Approach to Linear Model Analysis for Multivariable Turbine Engine Control Design 

C. A* Skira, Air Force Aero-Propulsion Laboratory, and R. L« DeHoff, Systems Control, Inc* (Vt.) 

SESSION 2. Theme Session A: Inverse Nyquist Array 

Chairman: B. Lehtinen, NASA Lewis Research Center 
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9:00 The Inverse Nyquist Array Method 

H* tU Rosenbrock and N* Munro, University of Manchester, England 

10:00 Insight into the Application of the Inverse Nyquist Array Method to Turbofan Engine Control 
H. A. Spang, III, General Electric Research and Development Center 

10:30 BREAK 


SESSION 3-1* Trans/fei* Functions I SESSION 3-2* Alternate Methods 

Chairman: S. Kahne, Case Western Reserve Univ* Chairman: J. Gibson, Texas A & M University 


11:00 Multivariable Design Problem Reduction to Scalar 
Design Problems 

B* D. 0, Anderson and N* T* Hung, University of 
Newcastle, Australia 

11:30 The Multivariable Nyquist Array: The Concept of 

Dominance Sharing 

G. G* Leiningcr, University of Toledo 

12:00 Input Compensation for Dominance of 
Turbofan Models 

R. H. Schafer and M* K. Sain, University of 
Notre Dame 


A New Frequency Method for Multivariable Systems 
R. DeSantis, Universite de Montreal 


Performance Analysis of Stochastic Linear Control 

Systems: A New Viewpoint 

S* R. Liberty, Texas Tech University 

An Automatic Depth and Pitch Control System for 
Submarines 

V. Nitsche, K. Luessow, and G. J. Thaler, Naval 
Post-Graduate School 
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12:30 LUNCH 
/ 

SESSION 4. Theme Session B: Complex Variable Methods 

Chairman: N. B. Nichols, Aerospace Corporation 


2:00 Complex Variable Methods for Multivariable Feedback Systems Analysis and Design 

A* G* J. MacFarlane, B. Kouvaritakis, and J. M. Edmunds, Cambridge University, England 

3:00 The Characteristic Frequency and Characteristic Gain Design Method for Multivariable Feedback Systems 
B* Kouvaritakis and J* M* Edmunds, Cambridge University, England 


3:30 BREAK 

SESSION 5-1* Transfer Functions II 
Chairman: B. Doolin, NASA Ames Research Center 

4:00 Linear Multivariable Control — A Problem of 
Specifications 

Z. V. Rckasius, Northwestern Uni*rnrsity 

4:30 Linear Multivariable Synthesis with Transfer 
Functions 

J* L. Peczkowski, Bendix Energy Controls Division 
and M. K. Sain, University of Notre Dame 


SESSION 5-2. Spectral Methods 

Chairman: R. DeSantis, Universite de Montreal 

Stability and Ilomotopy 

R. Sacks, Texas Tech University and R* DeCarlo, 
Purdue University 

A Compensation Procedure for the Desensitization 
tff Multivariable Regulator Eigenvalues 
P. J, M* Martin, R. K* Cavin, III and J* W. Dowse, 
Texas A M University 


r 


5:00 Application of Frequency Domain Multivariable 

Control Synthesis Techniques to an Illustrative 
Problem in Jet Engine Control 
L* G. Hofmann, G. L. Teper, and R. F* Wltftbock, 
Systems Technology, Inc* 


Linear Multivariable Synthesis by Eigenvalue/ 
Eigenvector Assignment 

S* Srinathkumar, NASA Langley Research Center and 
R* P. Rhoten, Oklahoma State University 


5:30 NO HOST COCKTAIL PARTY 


ORIGINAL PAGE IS 
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Friday, October 14 

SESSION 6. Theme Session C: Regulator Methods 

Chairman: S, Brodsky, Office of Naval Research 


Alternatives for Linear Multivariable Control 

N* Munra and 3. Hirbod, University of Manchester, England 

The Systematic Design of Linear Multivariable Control Systems for the Servomechanism Problem 
E. J. Davison and W* S. Casing, University of Toronto 

Linear Multivariable Control Design Based an Asymptotic Regulator Properties 
C. A- Harvey and G. Stein, Honeywell Systems and Research Center 


BREAK 

* 

SESSION 7-1. Modelling 
Chairman: K* Wozny, 

National Science Foundation 

A Direct Method for Obtaining 
Nonlinear Analytical Models of 
a Jet Engine 

R* J. Leake and J. G. Comiskey, 
University of Notre Darae 

A Multi-Time-Scale Design 
Approach for Jet Engine Control 
Systems 

A, .J* Calls e and B. Sridhar, 
Dynamics Research Corporation 

BLS as an Alternative to Linear 
Control Systems 

R. R. Mahler and V. R. Karanara, 
Oregon State University 


LUNCH 

SESSION 8-1* Output Feedback 
Chairman: E* M. Cliff 

Virginia Polytechnic Institute 

Output Feedback Regulator 
Design for Jet Engine Control 
Systems 

H, C* Merrill, NASA Lewis 
Research Center 

A Classical Root Locvs Design 
Method for Multivariable 
Systems in State Space Form 
G> K. Lee, University of 
Connecticut, M* Sohrwardy, 
Ruhr-Universitat, and D. 
Jordan, University of 
Connecticut 

Output Control via Matrix 
Generalized Inverse 
R. J, Miller, D* L* Powers, 
and V. Lovass-Nagy, 

Clarkson College 


SESSION 7-2* Model SoVla&ing 
Chairman: W. R* Perkins, 
University of Illinois 

Active Maneuver Load Control 
for a Control Configured 
Airplane 

N. C. Weingarten and E* G* 
Rynaski, Calspan Corporation 

A Parameter Optimization Method 
Applied to Engine Control 
System Design 

Y. Cheng, NASA Dryden Research 
Center 

Model Algorithmic Control 
A, Rault, J. Richalet and J. 
Papon, ADERSA/GERBIOS, France, 
and R* Mehra and W, C. Kessel, 
Scientific Systems 


SESSION 8-2* Additional 
Approaches 

Chairman: I. Rhodes, Wash. Univ. 

Optimal Open Loop Compensator 
Combined with Riccati Feed- 
back Compensator Control . * . 

R* Froriep, D. Joos, G* Kreis- 
selmeier, DFVLR, West Germany 

Observing Partial States for 
Systems with Unmeasurable 
Disturbances 

S, H, Wang and E* J. Davison 
University of Colorado and 
University of Toronto 


On the Design of Accurate 
Observers 

S. P* Bhattacharyya and 
I, G. Trindade, Universidade 
Federal da Rio de Janeiro 
and Universidade Federal 
Fluminense, Brazil 


SESSION 7-3* Comparisons 
Chairman: X. E, McDonald, 

Los Alamos Scientific Laboratory 

Quasi-Upper Triangular Decomposi- 
tion* Applied to the Linearized 
Control of a Turbo- fan Engine — 
Preliminary Results 
W. E. Holley, Oregon State 
University 

Reliability Considerations in 
Decentrally Controlled Multi- 
Variable Systems 

F. N. Bailey, E. B. Lee, and M. K. 
Sundareshan, University of Minnesota 

On Alternative Methodologies for 
the Design of Robust Linear 
Multivariable Regulators 
R. G* Kwatny and K* C. Kalnitsky, 
Drexel University and TASC 


SESSION 8-3. Design 
Chairman: E* C. Tacker 

University of Houston 

A Conceptual Design Approach 
Using Feedforward Plus Forward 
Compensation 

N* H* KcClararoch, University of 
Michigan 

^ Computer Aided Design of Control 
System via Optimization 
David Q* Mayne, Imperial College, 
London, England 

Design of Linear Multivariable 
Control Systems by CIP 
L, L* Gresham, J* R, Mitchell, 
and W# L. McDaniel, Jr* 
Mississippi State University 
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. "-.ABSTRACT 

. » Jhis work is *a comoarison of successively more 
^■comprehensive simulation models of an F— iQO-like 
: tTurhofan jet engine. A large and elaborate com- 
-..v^—puter program called DYNGEI* , developed over a 
,. ■-•jr^meber of years at NASA lewis Research Center, is 
“ 1 ^employed as the most comprehensive model for an- 
■. - .^alyzing steady-state and transient performance for 
* • -^control studies. This model employs many block 
Jafca maps and includes about 25 states. In order 
e£o perform optimal control studies, low order non- 
-linear analytical, and linear models have been 
/eloped. This paper reports on the details of 

* .these models and presents experimental data on 
.^tttheir relative performance. 

--^£HmiODUCTIOH' _ ' *. 

„u7.vTn .this paper we -consider the determination of a 
■\ -^Simplified nonlinear analytical ncdel for a two 
■spool turbofan jet engine. A large and elaborate 
(about 4000 FORTRAN statements) generalized engine 
. . . .^simulator called DYNGEN [1,2] ceded with represent— 
"dative block data maps, design parameters, and two 
^i-^spool operation is taken as the principal object 

• to he approximated, First we present the various 
- — ^models and then performance comparisons are made. 
,/-+5he models considered have been enumerated as fol- 

jlo vs« 

X 

i Model 0. Ih* actual jet engine (hypothetical.) 

Jfodel 1. The DYNGEtf simulator, coded with data pre- 
i sumed to have been taken from experimental measure- 
* —’Kents on Model 0« This model solves nore than 16 
nonlinear differential equations and uses data maps 
• JRnd thermodynamic tables which cannot be expressed 
.analytically, 

■ Hod el 2* This is an analytically expressed set of 
*5 nonlinear differential equations plus about 20 
’ * static equations expressing the relationship he- 
.t--£t?ecn various engine variables. The main task in 
; this project: was to determine this model, 

I 

Hod el 115, This is a normalized 5 th order linear 
*■ ••'*^®odel- which is obtained numerically frem Model l r _ 

’ — -Rising an experimental version of a program [3] 
^bcing developed at NASA lewis Research Center, 


rr^odel 113, This is a normalized 3rd order linear 
T.-^scodel* obtained by a hand calculated order reduction 
-«o£ Hod el 115. 

. -Model 215. This is a normalized 5th order linear 
-^anodel obtained by taking partial derivatives of the 
analytical Model 2. . 

Model 2L3. This is *& normalized 3rd order , linear 
■ -..model obtained* by a hand calculated order reduction 
of Model 215. 

. ^-ANALYTICAL MODEL 

■-In this section ve discuss and present Model 2., a 
k * jrfiimpllfied nonlinear analytical model of the jet 

! -» -engine. There are several reasons why it is desir- 
4 -nblc to Ixave such a model. First of all, one likes 
* io-sEe the basic nonlinear relationships between 
the engine variables in order to gain insight to 
■ - their dynamical and static behavior. Secondly, 

. such a model is invaluable in the application of 
■ .optimal control techniques [4] to engine control 
, system design. In the third place, if the model is 
I reasonably 1 accurate, it can. be employed as a fast 
r inexpensive nonlinear engine simulator for the 

evaluation of linear and nonlinear control strategies. 
Finally, linear models obtained by partial differ- 
entiation of this model tend to have more structure 
(zero entries in the A5CD matrices) than those ob- 
tained numerically, which gives the linear designer 
more insight* These linearizations then serve as a 
back tip to compare with the numerical linearisations . 


Jh the determination of Model 2. as an approximation 
of the DYNGEN Model 1. , theoretical relationships 
tdeveloped in [5] , [&]. and [7] were employed as a 
starting point. Certain simplifications suggested 
in IB] were used; and linear, least squares, ex- 
ponential, and polynomial fits to the Model 1. data 
About the chosen de3ign point were made* 

A 2 input* 5 state, 2 output model was developed 
with the following variable designations* 

*« fuel flow (WFS) 




nozzle area (Ag) 


! •OCjj'** - compressor rotor speed (K ) 



fun rotor speed (N ? ) . •• * 

burner exit pressure (P^) 
sfter burner exit pressure (P^) 
.high inlet energy (U^) 
thrust (PC) 

high turbine inlet temperature (T.) 


T n + 214.2732 CST - 48.0 (A d - 2.948255) 


JK^ESICN ^ I 21 /^ 2 10070 /T 21 /518.b68 


2,945255) 


3Che model is completely detemined by the .following 
'Specifications. 

. ^Constants 

-.3 ~ AJ - 773.26 M - 20.71175 

- 2.948255 £^.-.24 

_R - KA - -.0252 Cp p - .24 

~T* - 1.4 .-C^ - .20279 

T 2 - 1. Cgg T - .22589 

r-Xj - 318. L08 . ..Cp LT 1 .27938 

J_ - FttffiP - 3.8 4 - EEC -.16 - FCBLC 

C • VAC 

£1 — 'TMILP - 4.5 a. -"TCBLDU - .208 


6.20568 + .0129774 T, 


0185376 P 


3.516739 CNF - 63.916 


- 137.54 - 457.987 C!C + 564.325 CNC 
= 388.113 CNC 3 
, - 6.492 - 4.9747 CNC 
= 26.43184 - 69.0484 CNC ■+ 109.7234 
CNC 2 - 36.5756 CNC 3 


COMB 


CMAX 


-AFBN 

_JN_ -DESIGN - XNHPOS = 10070 


DESIGN - XNLFOS - 9651 


CMAXCKAX 


■DcbIkp Values 


EG - 13431.02 
T. « 2982.04 


N c - 11899.2 
K f - 9873.94 

- 23.9299 

- 586.467 

* - 2.55007 


(16) VA3 « (1- t>)KAC - .84 VAC 

(17) VC50 - 301.957 P4A' t "T7 — 


3.948255 


’(18) VG4 - VC50 - BsiVAC = VG50 - .116164 VAC 
^(19) VG55 « VC50 + ytVkC = VC50.+ .01056 VAC 


106.002 + . 86154T, 


10458CNC /T _./ 


55 414.532 7 

2.01365 
1121.784F 7 Ae 


+ c^^c55(t 50 - x 55 )] 

^ [T, (WA3 + WFB — WG, ) ] 


02951 WG7 /1934.415T. 

♦2116. 217A, 


(.5397SP--1) 


, [T, {VGA - WFB - WA3) 

COMB* A 4 

+ Y*(T.VA3 - T.VGA + T. (1 + n)VFB}] 


ypnllnear Functions Required for State 
end Outputs 


Figtrre 1 * Tvo-spool, two-*treaa 
turbofan engine 


H^JLSIGX 


dt 

v 1 

.COMB 

< P 7 

: . 

dt 

y 1 1 

AFBN 

"4 



< 

' so 4 

> < • 

1 

.k_i 

A_ 



a. 254 __ —.047955 
-5.9244 * .0048499 

97.919 -74.905 . 


Air flow 
Gas flow 
Pressure 
Temperature 
Internal Energy 


..•79723 
- —36146 
.... — 39*261 

0 1.461 

-.97908 —.01486 

. . 1.2138 

-15978 

•74.976, -7.136*1.1717, 


•-Olumbcra refer to the location within the engine 
t._(Hgure 1). All variable names correspond to those 
-tvhich are used in DYNGEN 


63298 


•llhc following linear no dels were obtained, 
V-vlth thrust, (FG), as the output if there is only 
^one output. 


JEIPERIMH3TAL COMPARISONS 

Figures 2 and 3 represent the nonlinear steady 
-atate operating lines of the fan and compressor 
respectively, as both nozzle area and fuel flow are 
changed. Ibis is a derailed nap in the vicinity of 
our nominal engine design configuration. The design 
point for model 2 differs from that of model 1 by 
-.less than 1/21. 


-1.152 1.448 
-1.991 1.071 

86.807 76.91 
-88.69 27.83 
—37.40 —246.7 


Model 2 
•vary nozzle 
area 


I —.8594 -.1397 .6672 1.167 .1236 

0555 91 .00656034 -.C018374 .0135393 .85391 

- f -.102766 1900938] 

^ l— .013839 .02085 6J 

••iZigenvelues: -251 ± 23J, -96, -5 + 0.8j 

~^Mod*»l 1L3 


Pressure 

Ratio* 

(*21**2?° 


Model 1 
^ry noz: 


Model 1 > 

vary fuel flow^x^y 


^riodel 2 ^ 

vary fuel flow 


-2.4307 -.70837 
3.8261 -4.9579 
2.4466 140.5 


-.81149 

-1.7235 

•94.982 


Z<4 2/0 2 2Z 

Fan Air Flow (WAT) (lbo/sec) 
Figure 2. Steady State Fan Map 


’^1543 .013382 1.333J °[.2346 .,87572) 

Eigenvalues: -92, -5 + 5j 
Model 2L5 

’-12.5487 -1.59279 3.58369 .237360 1.70112* 

.833048 -5.61346 1.64496 0.143805 1.30396 

lm 671-604 387.711 -392.675 -26.1602 160.52' 

-104.155 21.1351 64.9255 -67.8031 2.63144 

. 50.9527 -55.8546 -81.2047 -7 *.47450-105*743 

0.0 1.40762“ 

! 0.0 .758167 

1.28129 -122.314 

0.0 -48.9280 

149.210 -3.09196 


Model 2 

vary fuel flow 


r Model 2 
vary nozzle 

area 


/ Model 1 
vary fuel flow 


'Compressor Air Flou (VAC) (lbn/sec) 
Steady State Compressor Map 


r o.o 

0.0 

• 0.0 

1.46096 

o.ol 

fo.o 

0.0 

0.0 

0.0 

i.oJ 


-x 00259 

•355326 

. ;2H6 

-.316176 

12.54 

-13.7828 

—.6201 

— 99.388 

157.783 

6.84396 




the tine responses of various states due to 51 
v~*tep Inputs in fuel flow and nozzle area are shown 
In Figures 4 and 5 respectively. Only the nonlinear 
models are represented. Although the deviations of 
r the 6tates fron their design point differ by as 
•rauch as a factor of two, the actual values (design 
-plus deviation) of model 2 remain within 12 of the 
.itaodel 1 values. For strps in fuel of minus 202, 
-tthe-acaces remain within £2 of each other. 


Model 2 


Model 1 


.• Change 

.Off 

iDcoign 


Time (sec) 

• -*:F!igiire 6. Thrust Time response — 

» Fuel Step Input 

-A comparison of the linear system frequency 
.^responses of nozzle area to thrust, and fuel flow 
-to thrust, is shown in Figure 7. Models 1L5 and 
.2L5 match extremely well for the fuel input, but 
not so well for the nozzle area input. Even the 
Tower order model 1L3 is closer in this case. 


Model 1 


N Model 2 


~T iwt (sec) 

figure 4, State lime Responses-Fual Input 


V , Model 2 


.U, Model 1 


P, Model 2 


Frequency Responses— Inputs 
* to Thrust 


P- Model 1 


SUMMARY 

To the best of the authors' knowledge, no non- 
linear analytical dynamic models of a tvo-spooi, 
tvo-atrean jet engine have ever appeared in litera- 
ture. Tndecd, it is the value of the development of 
such a model which is the cost important considera- 
tion in the evaluation of our work, l.e., that a 
.good nonlinear analytical dynamic model will provide 
* flexibility and usefulness wnich is non-existent 
. *in present non-analytical Jet engine simulations. 

«A1 though, sooe significant discrepancies exist, our 
^aaodcl yields results vnich.are .accurate to within 12 


* Tine (sec) 

State Tice Rcspooses-h’ozzle 
Area Input Step 


The effects on thrust by a 52 step in fuel flow 
•-ere shown In Figure 6. Model 2L5 yields results 
N*vhich are quite close to model 2, although, both 
• Are significantly different from model 1. 




j.me* r the design point, and which degrade to an ac- 
curacy of approximately 8X with a drop in fuel flow 
of 20X4 The frequency response for fuel inputs of 
linearized models is also in close agreement. 

-In conclusion, ve are encouraged by our overall 
•-progress towards the development of the analytical 
*«Bodel, however, ve feel that more work is needed 
.--to Xurthcr improve the accuracy of our model. 
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Abstract 

Following Eosenbroc^s ideas regarding the advan- 
tages of dominance in linear multivariable control 
systems, a new graphical technique is used for Che 
design of compensators that achieve dominance. The 
technique is illustrated with an application to the 
problem q£ designing compensators for a linear tur- 
bofan-enginc model. The resulting design is put in- 
to perspective by examining it in the light of two 
other multivariable frequency-domain methods. One, 
MacTarlnnc: 1 s method of characteristic loci, is used 
to realise a final design for stability and low 
Interaction- The ocher is a direct technique based 
upon the algebraic expansion of the determinant of 
the return difference in terms of it's elcmcurs* 
Results from simulations carried out on the NASA 
DYMGKM software arc included. 


V" 1 ■ Introduction 

Recent years have witnessed a renewal of interest in 
frequency domain design methods for linear multi- 
variable control systems. The preponderance of 
these ideas are closely related to classical Nyquist 
constructions on the determinant of return differ- 
- ence. In this paper, we use three such methods to 
design a compensator for a two-input, five-state, 

. two-output linear model of a modern two-spool turbo- 
fan jet engine obtained from the DYNGEM digital jet 
engine simulation* 

Rosenbrock [1] has related the classical Ny^uist 
construction on the determinant of return differ- 
ence Lo corresponding classical constructions on the 
diagonal elements of the return difference — pro- 
vided these diagonal elements "dominate 11 their rows 
or columns in an appropriate manner* Focussing the 
design interest on achieving dominance in this 
sense. Section 3 presents a new graphical technique 
to help with this aspect of design* Next, Section 
4 utilizes the general! ted liyquist plots to obt iin 
an acceptable compensator design. The ideas of 
generalized Nyquist plots were introduced by 
MacFarlauC [2], who related the determinant of the 
• return difference to it T s spectrum when regarded as 
*n appropriate linear operator. 

In Section 5, we utilise a direct technique which 
emphasises the algebraic relationship between the 

Hi is work, was supported in part by the Rational 
Science foundation under Grant ESC 75-22322 and in 
part by the National Aeronautics and Space 
Administration under Get ant: MSG 3048* 


dements of the return difference and it’s deter- 
minant* Typically, when it achieves user satis- 
faction, this method does so with greater speed, 
and fewer concepts, than it f s competitors* The jet 
engine model is introduced in Section 2, which also 
establishes the notation for succeeding sections* 
Finally, in Section 6, we give results of simula- 
tions to evaluate the performance of the system* 


2. Jet Engine Model and Return Difference 
Determinant 

The linear model used far the study is based upon 
data obtained from a DYEGEN simulation. Xt is 
specified by the equations 

x = .Ax + Ru (1) 

y * Ox + Du (2) 

where x, u, y denote the state, input and output 
Vectors respectively. The inputs are fuel flow 
and nozzle area; the five states are compressor 
rotor speed, fan rotor speed, burner exit pressure, 
afterburner exit pressure and high pressure turbine 
inlet energy; while thrust and high pressure tur- 
bine inlet temperature constitute the two outputs* 

Ke next consider the problem of designing cox: tr oi- 
lers for the plant. The underlying feedback con- 
trol scheme is shown in Fig, 1. G(s), the plant. 



Fig, 1 Feedback Control Scheme 


represents the jet engine model, that is, 

G(s) *= C(sl - A)~^B + D. (3) 
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Xj and y lie inside and dashed circles if the 
acceptable region Is outside. 

If the above procedure in repeated over a range of 
frequencies for each vow of the system, and the 
circles of intersection drawn, a plot describing 
the acceptable values of the complex number x^ 4- 
Jy^ for each frequency results. In this way, the 
acceptable range of the function S.(s) such iliac 
the i th row of K(s)C(s) is dominant is described* 

Thc^analysir. of the CARDIAD plot for a given row 
of C(s) proceeJs as follows. If the origin of the 
plot is contained inside all solid circles and is 
excluded by all dashed circles, the row of C(s) is 
dominant uncompensated* If the row of G(s) is not 
dominant uncompensated , the CARDIAD plot is next 
checked to see if there is a constant entry S 4 that 
will make K(s)G(s) dominant at all frequencies. 

For this to be the case, there must be a point on 
the real axis that is included in all sella circles 
and excluded by all dashed circles. 

If there exists no constant £. such that the i 1 ^ 
row of K(s)G(s) is dominant al all frequencies, the 
CARDIAD plot is used as a guide to design a fre- 
quency dependent S.(s) that will achieve dominance. 
This is accomplished by realizing a function 3^ (r) 

whose value at iw lies inside the circle associ- 
c 

ated with the san? frequency in the CARDIAD plot if 
that circle is solid, or outside if that circle is 
dashed. This approach is illustrated by consider- 
ing the DYNGEN problem. ; 


Central to the application of Nyqulst type Ideas to 
multivariable systems Jn the return difference ma- 
trix, which In this case becomes [ I !G(s)K(s) ] . 

It's principal use arises from the relation of the 
closed loop characteristic polynomial (CI.CP) to the 
open loop characteristic polynomial (01.CP) which 
can be stated, In the manner of [1], as 


where equality is understood up- to a real constant. 

Of primary concern here is the behavior of dct(HGK) 
for values of s on the standard Nyquist contour 
(SMC), which encircles the open right half plane 
clockwise with indentations into the left half plane 
around poles and zeros on the imaginary axis. In 
practice., plols arc made for values of s on the 
positive imaginary axis. Stability can then be de- 
termined from plots of dct(I+CK) in conjunction with 
knowledge of the open loop characteristic polynomial. 
Also interesting, of course, is the use of such plots 
to aid in the choice of a suitable K(s). 


3. CARDIAD Plots and Dominance [3,4] 

The CARDIAD (Compensator Acceptability Region for 
DIAgonal Dominance) plot is a graphical approach to 
the problem of choosing a compensator that will 
achieve system dominance. A system is said to be 
row (column) dominant [1] if the magnitude of each 
diagonal element of the open loop transfer function 
matrix is greater than the sum of the magnitudes of 
the off diagonal elements of the row (column) at all 
frequencies. In the 2x2 case being considered in 
this paper, the dominance condition reduces to the 
magnitude of the diagonal element being greater 
than the magnitude of the off diagonal element of 
the row (column) . Consonant with the Rosenbrock 
approach, the CARDIAD plot analysis is applied to 
the inverse of the plant G(s). As a notational 
point, the inverse plant transfer function matrix 
will be donotedby G(s) and the inverse pre- 
compensator by K(s). ; 

The specific application to the jet engine design 
problem involves trying to find a compensator K(s) 
such that K(s)G(s) is row dominant* Without loss 
of generality, the form of K(s) will be restricted 

to 


mm 


where 


If G(s) and R(s) are each evaluated at a frequency 
J[w i^the equation for dominance of the i r ^ row of 
K(s)C(s) becomes a function f . (x^ , y )vhich describes 
a paraboloid in three-space. The intersection of 
this paraboloid and the complex plane is a circle 
which is the locus of the values of and y. such 
that the magnitude of the diagonal element of the 
ith row* of £(s )G(s) is equal to the magnitude of 
the off diagonal element of the row. Minima and 
iiaxlna analysis of the function f. revealr that 
values of Xj and y^ on one side oi the circle will 
inake the. system dominant, whereas values which lie 
on the other side ot the circle will not. In the 
CARDIAD plots, this differentiation is made by 
drawing u solid circle if the acceptable values of 


Fig. 2 CARDIAD Plot Row 1 
Uncompensated 


The initial CARDIAD plots of G(s) indicate that 
row 2 of C(s) is dominant uncompensated since the 
plot consists only of dashed circles, which all 
exclude the origin. The plot for row 1, however, 
shows that this row is not dominant uncompensated 
and also that there is no constant entry in the off 
diagonal clement of row 1 of K(s) that will make 
the row dominant at all frequencies. This is 
easily seen since all the circles Jn this plot are 
colid and there is no point on the x axi-s that is 
included in all the circles. Moreover, the plot 
hints that there will be difficulty finding a 3,(s) 





Fig. 3 CARDIAD riot Row 2 
Uncompensated 

that will make this row of £(s)G(s) dominant because 
of the complexity of the plot and the snail radii 
of the low frequency circles which necessitate a 
very close fit. 

To facilitate the process of finding a compensator 
that will cake K(s)G(s) dominant, the system was 
first precompensated with 








Space limitations do not allow the CARDIAD plots of 
K-G(s) to be included, but the new ^lots arc the 
same shape as the CARDIAD plots of G(s) with two 
major changes. The row 1 plot of K G(s) is the 
same shape as the row 2 plot of G(sj with dashed 
circles changed to solid circles. Similarly, the 
row 2 plot of K.G(s) is the same shape as the row 1 
plot of C(s) with the solid circles changed to 
dashed circles. 

The. problem of finding a K^(s) such that ic o (s)K 1 G(s) 
is^dominant is now simplified. Since tow 2 of 
K^G(s) is now dominant uncompensated, the off diago- 
nal term in the second row of ^(s) is left a 2 ero, 
with the provision that if it later proves helpful 
in compensation, the entry may be chosen to be any 
constant that lies outside all of the circles. To 
make row 1 of R 1 6(s) dominant, the off diagonal 
entry in row 1 of r.ust follow the semicircular 

path through the complex plane described by the 
CARDIAD plot for this row. A fit was made to this 
shape and the resulting K^Cs) was 

9.4798 + 0. ?494s 


* 2 (s) 


1.2339s 

1 


(B) 


The CARDIAD plots of K ? (s)KjC(s) arc considerably 
*»rc complex than the previous plots. The plot for 
row 2 shows that the row is dominant at all fre- 
quencies since the origin is included by all solid 
circles anJ excluded by all dashed circles. The 
CARDIAD plot for row 1 shows that the row is clearly 
not dominant at all frequencies. Dominance is lost 
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Fig. 5 


Row 2 


CARDIAD Plot 
Compensated 

at v«10, is regained as v=90,« and is lost again at 
v®700. It is perhaps possible to.find a better 
choice of J^Cs) that will make row 1 dominant at 
all frequencies, but the dominance achieved by tin* 
above K^(s) proved to be sufficient. 

An interesting feature of the CARDIAD plot is 
Illustrated in the final plot for row 2. Close 
analysis o' this plot shows that there are three 
occurrences of solid circles changing to dashed 
circles or dashed circles changing to solid. When 
these transitions occur, the paraboloid is invert- 
ing and the circle of intersection degenerates to 
a line. These lines occur when the other row 
changes from being dominant to not dominant or vice 
versa. Thus, each change in dominance of row 1 
causes a change in the type of circle being drawn 
In the plot of row 2. , 

I 

i 
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4. Design Using Characteristic loci 

Another approach to design, due to A.C. J.MncFarlane, 
uses the locus of the eigenvalues of G(s)K(s), call- 
ed the characteristic loci (C.L.), for values of s 
on the SNC. This method is based on the relation 
of the determinant of the return difference to 
eigenvalues of C(:*.)K(s). In order to assess sta- 
bility from the C.L. plots, for s a jv, v positive, 
one must count the clockwise encirclements of the 
critical point (-1,0) inade by the C.L. plots and 
sun all these up. The closed loop system Is stable 
if this sum equals -po, .here p--> is the number of 
zeros of OLCP enclosed by the SN'C. As an approxi- 
mate measure of interaction, vc compare the eigen- 
value plots with plots of the diagonal elements of 
\ Q(s)=G(s)K(s) . For a nouintcracting system with 
Q(s) a diagonal matrix, these would be identical. 

In our design example, Q(s) is a 2x2 matrix, and 
therefore we will be looking at plots of two cigen- 
• values A^(s), A^s) anc * tvo diagonal elements 
q 22 (s). 

First, an examination of the C.L. plots of the un- 
compensated system revealed that, without compen- 
sation, the closed loop system is unstable. The 
plots arc not included due to lack of space, but 
conclusions drawn from them are given. Control pro- 
blems for the uncompensated model were complicated 
by the existence of considerable interaction, and 
large gains at high frequencies. An additional 
difficulty was that one of the eigenvalues was nega- 
tive at zero frequency. This tended to limit the 
response speed the closed loop system. It ap- 
peared on the C.L. plots that, from a stability 
vicwpojnt, the frequency range of interest is in the 
* vicinity of 10 rps. This gives justification for 
use of the compensator given in the previous section. 
As a practical nutter, our goal is to achieve as 
rapid a response as possible to a step input, with- 
out suffering any overshoot. Heavy emphasis is 
placed also on steady state accuracy. 

To remove the right half plane pole in K ? , we choose 
arbitrarily as diag (1/s, (-l+1.2359s)/s) . Hie 
resulting becomes 


r(s> 


[i 


-1+1. 2359s 
9.479S+0. 2494 


] 


(9) 


The diagonal nature of K^(s) docs not affect domi- 
nance. Moreover, an examination of the (1,1) and 
Ot 2) elements of C(s) reveals that if the 0 in K(s) 
is changed to 9, vc can significantly reduce the 
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high frequency magnitude of q^(i;) while simultane- 
ously boosting the low frequency magnitude. This 
is in accordance with the freedom specified for L 
previously. The modified K(r.) gives rise to the 
plots of Fig. 6 and 7. 

Since dominance is not affected by diagonal compen- 
sators, the problem becomes that of independently 
shaping q. ^ and q^ by means of single loop tech- 
niques. In order to reduce high frequency gain in 
q^p without appreciably affecting low frequency be- 
havior, we use lag compensation. A little bit of 
cut and try led finally to K, (s)*diag. (0.44E-4 , 
(-s-0. 1 ) /(2000s+10) ) . The plots corresponding to 
lWC^K^K^K, are not shown. The and q plot 
is essentially that of A in Fig. 6 scaled by a 
factor of 0.44E-4. Similarly, the plot in- 
verted and scaled by 0.0005. (By inversion vc 
mean reflection through both axes.) The plot for 
*2 Is shown in Fig. 8. 

5, The Direct Method of Analysis 

Direct methods of multivariable Nyquist analysis 
concern themselves with the algebraic relation- 
ship between the elements of return difference and 
its determinant. For an NxM return difference, 
the most basic of these relationships is 

det(I+CK) » 1+1 {£ ixi principal minors of CK}. 

1=1 ’ (10) 

For the example of this paper, (10) takes the form 

l+{ ( G 11 K 11 + G 12 K 21 )+(G 21 K 12 + G 22 K 22 ) ) + dct(CK) . (11) 

In (10) and (11) we note tie advantage of minute 
detail and the disadvantage of nonrccursive con- 

struction. Considerable i* terest attaches to the 
removal of this disadvantage, which can be accom- 
plished by methods drawn from the results of exte- 
rior algebra [5], Consider the recursion (where 
tr denotes trace) 

°0 * 1 (12a) 

a r “ " r i » tr(GK) r-p < 12 b) 

T p=0 p 

for l<r<N. It can be shown that 
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Fig. 6 First C.L. plot for Fig. 7 Second C.L. plot 


CK 1 K 2 K 3 . 


for Gk^k^k^. 










dct(HCK) 


The direct approaches differ appreciably from meth- 
ods described in preceding sections, in that they 
address themselves directly to the image of 
|dct(lK:K) on the SN C, without any particular con- 
cern for such issues as dominance or interaction. 
Alternate insights accrue from such plots, which 
we ill ustrate here for the engine design example. 
All plots are drawn for the final return differ- 
ence as developed in the section:; preceding. 

Fig. 9 indicates the five constituents of a 
dct(HCK) plot as developed in (11), while Fig. 10 
presents the corresponding two constituents accord- 
ing to (13). Fig. 11 contains the total Nyquist 
plot, which is obtained by adding the individual 
curves in cither of the two prior figures. Re- 
vealed in this plot. Fig. 11 is a feature not so 
readily noticeable in the earlier plots, namely con- 
ditional stability. It appears, therefore, that the 
availability of a variety of graphical tools is in 
the multivariable case every bit as valuable as in 
the more classical, one-input, one-output situation. 

It is readily seen in Fig. 11 that the plot encir- 
cles -1-fjO twice in a counterclockwise direction. 
Therefore, the system is shown to be stable because 
the open loop characteristic polynomial has a double 
rcro at the origin. 

Further exploratory studies of direct methods as 
design aids arc available elsewhere [6). 


Fig. 10 The Nyquist 
Plot of the Elements 
of Det(I+CK) Accord- 
ing to Expansion, Eq. 

(13)- - 


Fig. 9 The Nyquist Plot of 
the Elements of Det(I+CK), 
According to Expansion Eq. 
( 11 ). 


6. Simulation Results 

Closed-loop time responses were obtained both by 
using the linear model simulation and by implement- 
ing the compensator on DYNGEN, a jet engine simula- 
tion program developed by the NASA Levis Research 

Center [7]. 

DYNGEN is a versatile digital program which ana- 
lyzes steady-state and transient performance of 
turbojet and turbofan engines. It uses a sixteenth 
order system to model this example, and solves 
the state differential equations by a modified Euler 
ncthod. The user need only supply appropriate com- 
ponent performance maps and design-point informa- 
tion, and then write the control subroutines. Im- 
plementation of the compensator required first order 
functions to perform integration and lead-lag com- 
pensation. 

The linear model used in this study was also obtain- 
ed from DYNGEN. By utilizing a special control sub- 
routine written by NASA, calle.d DYCABCD [8], models 
can be derived using whatever states the user de- 
sires. DYNGEN thus possesses the capability to de- 
termine linear models for the engine with any order 
up to sixteen. 

Fig. 12 shows the response of the linear model to a 
step input in the first channel. Thrust has a rise 
tim (10Z-90X) of 1.04 second:; with no overshoot 
occurring. High pressure turbine inlet temperature 
increases to a maximum of 0.105 at approximately 
0.9 seconds, Chen gradually decreases. 

Similar, and even better, results occur when the 
compensator is employed in the DYNGEN simulation 
using a one percent step. Thrust rise time is 0.88 


Fig. 11 Nyquist Plot of Det(I+GK) 


seconds, and the turbine temperature reaches a 
maximum increase of 0.097. • 

The linear and nonlinear responses were not in such 
close agreement for a step input in the second 
channel. The linear model shows turbine tempera- 
ture slowly ramping up (Fig. 13) as the change in 
thrust is held to a minimum. DYNGEN produces simi- 
lar results for the turbine temperature response; 
however thrust experiences a strong decrease bc^or? 
rising to zero. At this writing, it is believed \ 
that the five states chosen for the DYCABCD model 
dc not adequately describe local engine behavior 
for the second channel equipped with the present 
controller. 
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Fig. 12 Response to Seep in First Channel 
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Fig. 13 Response to Step in Second Channel 
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7* Conclusions ; 

This paper has demonstrated the usefulness of the 
new CARD I AD plot approach to designing compensators 
for complex plants. The DYNGEN simulation for a 
step in channel 1 has shown that acceptable re- 
sponses can be obtained using linear compensators. 

An ordered collection of these may make global con- j 

trol feasible. For steps in channel 2, conclusive 

evidence was not obtained. We suspect that this is ; 

due .to inadequacy of the linear model in describing 

the plant. This important factor of selecting an • 

appropriate linear model is often overlooked. But, 
as ve have seen, it turns out to be crucial in 
practical applications. 

The method of CARDIAD plots can be generalised to 
plants with more than two inputs and outputs by con- 
sidering a family of compensators with l f s on the 
diagonal and only one non-zero off-diagonnl term. 

As stated in [1), except for changes in the order- 
ing of Inputs or outputs, such a study is exhaus- 
tive. 
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SUMMA RY 

This paper describes a complete case study of the application of the theory of minimal de- 
sign to multivariable control of jet engines. The minimal design problem is approached from the 
viewpoint of polynomial modules, and computational experience with PL/I and FORMAC-PL/I software 
is discussed. The complete minimal design solution exhibits flexibilities not apparent in early 
industry studies, and a new approach to pole assignment can be used to advantage in this situation 

control system design in tne near term. 


1. INTRODUCTION 


One way to approach the design of linear multivari- 
able control systems is to express system specifications 
In terms of a desired closed loop transfer function ma- 
trix. A question which is often raised about such an 
approach is the practicality of making such a specifica- 
tion. Another, related, question concerns the possibili- 
ty of determining the existence of realizable compensa- 
tors to achieve the specification. When such compensa- 
tors do exist, there are the very practical issues cf 
giving a finite enumeration of them, of determining 
whether they have fixed poles, and of assigning one or 
more of the non-fixed poles. Of special interest, as it 
turns out, for the issue of pole assignment is the idea 
of minimality, in the state-space sense, of a proposed 
solution in the context of all possible solutions. 

This paper provides a thorough case stud} of such a 
design approach when applied lo realistic numerical mod- 
els associated with an F-100-like turbofan engine. 

Specifications are accomplished by means of the methods 
of linear optinal control theory, according to proce- 
dures already worked out in the jet engine industry, 
pie remaining tasks are addressed by regarding the design 
As a problem in free polynomial modules. A special fea- 
ture of the application lies in its attention to compen- 
sators of simple structure, with a view to the use of a 

t raded collection of them for the purpose of global 
ngine control. 

Section 2 describes the basic design problem, once 
pecif ications are made. Section 3 provides the dis- 

f ' ussion of the jet engine application, with particular 
ttention paid to the manner of making the specifications 
And to the formulation of the main design problem for 
the jet engine application. Section 4 explains how to 
gast the design problem in terms of free polynomial iod- 
ijles, and Section 5 describes floating point computation- 
al experience gained in applying extended precision PL/I 
software to solve the jet engine problem in the free mod- 
ule context. Section 6 outlines the corresponding 
experience associated with an exact rational calculation 
made with FORMAC-PL/1 software. 

The results of Section 6 show that considerably 
greater compensator design freedom is available than had 
been apparent from early industry studies. Using these 
results, n new pole placement design procedure based on 
Alternating multilinear algebra achieves in Section 7 a 
minimal pole placement solution not possible by those 
earlier industry methods. 

I I 

Section 8 closes with remarks designed to place the 
work in historical perspective, to reference the litera- -Beyond the basic LDP, several additional issues are 

ture, and to assess the aerits of polynomial methods for of practical importance. Among these should be included 

* This work was supported in part by the National Science Foundation under Grants CK-37285 and 
ENG 75-22322 and in part by the National Aeronautics and Space Administration under Grant NSG 3048. 


2. T1IE MINIMAL DESIGN PROBLEM 


Suppose that F is a given field. For the jet 
engine control problem, F is taken to be R, the field 
of real numbers; however, a great deal of the algo- 
rithmic nature of the discussion is more general than 
that, and is so stated. The set of polynomials which 
are of interest is F[s], namely those polynomials in 
the variable s with coefficients in the field F. The 
fact that F[s] is a principal ideal domain ring is well 
known, as is the equally pertinent fact that F[s] has a 
quotient field F(s). More intuitively, F(s) is often 
described as the field of rational functions in s hav- 
ing coefficients in F. ; 

The design problems of Interest in the sequel are 
conventionally stated in terms of F(s); however. 

Section 4 explains how such problems may be re-converted 
back to a corresponding F[s] form. 

Principal interest centers upon the minimal design 
problem (MDP), which can be described as follows. 

Let G: be a linear operator for finite-dimensional 

F(s)-vector spaces V and V 9 . G is regarded as 
realizable if its matrix is'proper. Now let C 7 : 
and G be given linear operators, where v., is 
also a finite-dimensional F(s)-vector space. MD^ con- 
sists in determining whether there are realizable linear 
roperators G which make the diagram in fig. 1 commute 
and, if so, to find one whose minimal realization is of 
least dimension among all such realizable operators. 
Intuitively, 
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(3) | 



;(i) a finite enumeration of all possible solutions, (ill 
jdc termination of any fixed poles In the matrix of G, and 
j(lii) methods for assigning the poles of C which arc not 
fixed. The answers to these questions resolve such 
{issues as the availability of solutions with varying 
•degrees of Integration. From a conceptual viewpoint, 
{these ideas are developed further in Sections 4 and 7, 
whereas the computational issues are discussed in 
Sections 5 and 6, | 

Next in order of presentation, however, is the 
statement of a minimal design problem for jet engine 
control. 

3. JET ENGINE APPLICATION 

In this section, we demonstrate the practicality of 
the minimal design approach in the context of jcr. engine 
{control. The basic plant is a version of the F-100 
jturbofan engine. Inputs are jet exhaust area and main 
jburner fuel flow; states are fan inleL temperature, main 
burner pressure, fan speed, high compressor speed, and 
afterburner pressure; and outputs are thrust ana high- 


J - j | r'y T Q6y + 6u T Rii + 6y T S6y)du 
0 

where superscript T denotes matrix transposition. The 
weighting matrices Q,R, and S arc listed in table 2. 

Table 2 

Weighting Matrices 

with Optimal Integral Control Solution 



I 



turbine inlet temperature. The linearized model approxi- 
mates the small signal behavior of these engine variables 
jin a neighborhood of 47° Power Lever Angle (PLA) . This 
corresponds to a point approximately midway between 

The plant 
in (1) and 


At this pMnt, a minimal design problem can be 
brought into play. The control scheme of fig. 3 is 
seen to be more desirable because it incorporates out- 
put feedback and enjoys the concomitant advantage of 
zero steady state error, even in the presence of plant 
parameter variation. 


is specified by th« 


Table 1 lists these matrices for our example. The 
attempt to design simple compensators for linear control 
pver a specified region 


strategy for global 
control of the engine using a graded collection of these 

Tabic 1 j 

State Description Matrices for Jet Engine (1LA®47°) 
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Fig. 3. Optimal Integral Control 

One relates the performance of the two control schemes 
by equating, in both, the Laplace transform of the 
variable 6u, as written in terms of the respective state 
variables. This leads to the following equations, which 
may be solved for L and H, the values of which have been 
listed in table 2. | 


Hint this is nothing but a form of the minimal 
design problem can be seen oy evaluating the 2x2 closed 
loop transfer function matrices T(s) and T r (s) for the 
two systems in figs. 2 and 3. In fig. 2 t 

| T(s) - P 1 (s) (si - r^s)) -1 M (6) 

where 

P.(s) - C (sI-A )“* B 4- D (7) 


We nc.*t examine how engine control specifications 
an Lc obtained from linear optimal control theory. In 
ig. 2, the compensators. 


Plant 


where 


Now, rewrite (6) using (4) as 
T(s) - -P, (s) (sI-P,(s)J 


Fig. 2. Linear Optimal Control 

nin matrices Gj and C^ are chosen with the 
inimizing the performance index of (3), 


The relationship between (9) and (11) now depends upon 
that between (R) and (10). 
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Comparison of P.j(s) and P 0 (s), with the aid of (5) 


and (7), establishes the equality of the tvo transfer 
‘functions T(s) and T'(s). 


We can then pose questions regarding the existence 
pf compensators other than H and L to achieve the same 
performance as attained in fig. 3, and wh.'t, if any, 
advantages such compensators would have over that scheme 
To do this, we consider tig. 4, which is a more general 
scheme of control based on fig. 3. 
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Wo remark here on a special feature of (20), in 
that it implies a common set of dynamics for both com- 
pensators. Suppose that a suitable solution of (22) 
leads us to matrices A , B , C , and B^, C^, and 


such that 


v-l. 


I l 


C(a) - C g (sl-A g ) B g+ D g 
K(s) - C k (sI-A v r 1 B k + D k . 


(23) j 

(24) | 


6z 





Plant 

(T\ - 

C(s) 

+/Q. 6u, 

(A * B , 
P P 

v y — 



C , D ) 

1 
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P P 


K(s) 
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In a typical situation, D * 0; and it can be shown that 
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kr-*y 
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Fig. 5 then shews how such a control scheme can be rea- 
lized. ' 


6z 






Fig. 4. Generalized Compensation Scheme 

pur objective is to design compensators G(s) and K(s) to 
achieve exactly the transfer function T(s) ■ T*(s) 
between 6z and 6y. This means that we must have 

(I + Ms))" 1 P, (s) » T(s) 


4 W/ ‘4 

Where we have introduced 


( 12 ) 


ena 


P 4 (8> 


P x (s) 


(I+K(s) Pj(s)) 


-1 




C (sI-A ) 
g Z 



Plant 

<v V 


c. D ) 
P P 


5x 


&y 


C(s) 


P 5 (s) 


( S I-A p ) 


(13) 
(14 ) ! 


From (12), we obtain the equivalent condition 


P 4 (s> 


(I+P^(s)) T(s), 


I 


(15) 


yhlch c.'a be restated in the manner 


Pjd+KPj) - ^ 


{x+p 1 (i+k? 5 ) 


-1 


G)T. 


(16) 


Fig. 5. Realization of Compensation Scheme 

Finally, we make the observation that (22) can be 
completely solved, and a minimal solution computed by 
algorithms given in the Appendix. The next section 
deals with the theoretical foundation of these algori- 
thms; and subsequent sections describe their application 
to (22). j 

4. FREE MODULAR APPROACH TO MDP 
In matrix form, the minimal design problem of fig. 1 


reduces to solving an equation 


T rom table 1, D is clearly invertible; and so the linear 
dynamical systeG has a unique linear dynamical inverse 
gystem P uhlrh up Hpcipnafn P Thuc MAI 


G 1 (s)G(s) 


G 2 (s) 


I 


(25) | 


which in turn can be written 


borne simplification can be achieved at this point if we 
take advantage of the fact that the matrix T*(s) has an 


[inverse T^s). Then (18) can be cleared in its right 


11 

C 

[ 

II 

r. 

o 


member so that 


which we designate P^. Thus (16) is 


Equivalent to 


G(I-T) - K PjPjT 


PjT 


(17) 


for the various realiza^e G(s), where G^(s) and G 0 (s) 
arc given. Section 3 provided a nontrivial illustration 


pf (25) in (19), where 


[(1-T T ) 


-T T pJr£] 


K 


T*T 

TP 
1 P r 


G 1 (s) 


[Pj (s) (T T ( s)-1) 


- rj«] 


(18) 


C 2 (s) 


I, 


(26) 

(27) 


and where the field F was R, the real numbers. 


The free modular approach to MDP is based upon the 
recognition that, as a set, 


F(s) e F[s J x F[s] , 


(P P (T T -I) 


T 

- P 5 ] 


(i9); 


f ow compare (19) with fig. 1, from which it becomes clear 
hat Qj is the identity map, or that our control system 

E inimal design problem turns out to be a version cf the 
inimal inverse system problem. Writing 
T 

- N(s)D -1 (s) , (20)| 

knd ,K (s) 


which, in turn, suggests that it may be possible to ex- 
press (25) in terms of F[s). A convenient way to bring 
this about, as illustrated in (20), is to write 


G(s) 


N(s) D -i (s), 


(28) ; 


lP*(s)(i T (s)-I) 


* P 5 (S>1= dU)* 


M(s) 


( 21 ), 


where N(s) and D(s) have their elements in Ffs). It is 
easy to see that every G(s) has representation in the 
form (28). Similar representations could be adopted for 
Cis), i * 1,2, but the presentation can be simplified 

it the choice ^ \ \ 

I M (s) 

1 i - l t 2 (29) | 


G jL (s) 


d^s) 


where N(s), D(s), and M(s) are matrices over R[s], and 
yhere d(s) e R[sj, we can put (21) in the polynomial form 


is made, with d,(s) c F[s], i * 1,2 and M^s) having ele- 
1,2. Equation (25) is then clearly 


fM(s) I — d (s) 1 ] 


N(s) 

P(s) 


( 22 ) 


mcnls in F(s], i 
the same as 


Equation (22) is a polynomial "kernel" problem, 
Equivalent to the design problem of fig. 4. By compar- 
ing figs. 3 and 4, we can trivially establish that a sol- 
ution to the problem does indeed exist. Our goal in the 
fcequel is to give a finite enumeration of all possible 
Solutions and to study their pole assignment possibili- 
tles relative to the structure of fig. 4. . J 

47 . 


M. (s) . H,(s) 

N(s)D“ A (s) - 


d^s) 


d 2 (s)* 


(30) 


which, in turn, is equivalent to 


[d 2 (s)M 1 (s) 


-d 1 (s)M 2 (s)) 


an equation written over F[s]. 


[N(s)l 

lo’cVlJ 


0, (31) 

For the jet engine 
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[problem, (22) corresponds to (31). Let 

(32) 


denote the 1th column of 
N(s)l 


Tlien 


4(6) - 


n t (s) 


d^s) 


D(s) 


[d 2 (s)M JL (s): -d 1 (s)M 2 (s)] t A (s) 


0, (33) 

and every candidate tc construct a solution (28) can be 
traced to such ^(s). 

Thus MDP is quite closely related to the homogeneous 
equation 

[d 2 (s)M 1 (s): -d 1 (s)M 0 (s)] t(s) « 0. (34) 

|The purpose of this section is to explain briefly an 
[appropriate algebraic interpretation of (34). This in- 
terpretation is based upon generalising the notion of the 
[n-dioensional F(s)-vector space F(s) n to that of a 
irank-n F[s]-module F[s] n . As a vector space, F(s) n 
•satisfies the usual axioms, with scalars taken from the 
>f ield F(s). As a module, F[s] n satisfies exactly Lhe 
came set of axioms, but with scalars taken from the 
principal ideal domain rin g F ( s ] - Despite this close 
(similarity, F[s]-modules do not behave in exactly the 
same way as vector spaces. But there is a class of them, 
known as finite-rank free modules , which have a great 
similarity to finite-dimensional vector spaces in that 
jthey have a basis , which can be defined in the usual way 
bsing concepts of span and independence. F[s] n , for 
[example, is said to be free on the basis 

{(0, 0, 1, 0, .... 0); i «= 1,2, . , . ,n) . (35) 

/ 

ith position 

Morph ism s of Ff s]-modules are defined analogously to 
jilnear operators on vector spaces; and, when domain and 
codomain are finite-rank free modules, the basis concept 
jl6 used in the usual way to define a matrix for the 
morphism. This, then, is Lhe interpretation to be given 
to the p x q matrix I 


[d 2 (s)M 1 (s): -d 1 (s)M 2 (*)l (36) 

in (34), namely the interpretation of a morphism 

M: F[s] <? -*• F[s] P (37) 

>f finite-rank free F(s]-modules. As a submodule of the 
finite-rank free module F[s]^ over the principal ideal 
omain F[s], the kernel of M is also free, and thus the 
olution to (34) is tantamount to finding a basis for 
his kernel. The process for calculating such a basis 
e provided by Algorithm 1 in the Appendix. 


If a basis 
t 1 (s). 


t 2 (s). 


t^(s) 


t t (s) 


where t 




c F 1 * 


. r 

J-0 




(39) 


4 0, and i * 1,2,..., A. Then the 


basis (38) is said to be reduced if the matrix 


lt l k : C 2 k 1 
1 * * *2 * 


** k J 


(40) 


has rank I, Algorithm 2 In the Apm-ndix fjlllCM a b.i.- is. 


* ) 
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With those notions, the MDP Algorithm In the 
Appendix solves MDP. The issue of pole assignment is 
taken up in Section 7. 

5. FLOATING POINT EXPERT KECE 

In view of the material presented in the previous 
[section, we are ready to take a closer look at (22). 

,Tiie matrix [M(s)I -d(s)I] turns out to be 2 . 2 x 9 matrix 
y>f polynomials in R[sj. Lack of space prevents us from 
reproducing all the numbers here, but fig. 6 shows that 
the typical element is a thirteenth degree polynomial. 

We also note the large variation in the magnitudes of 
the coefficients of the polynomials. 




] 


os); 

or Ker At has been computed, MDP solution then depends 
ppon a determination of whether these basis elements can 
be used, through (32) and (28), to construct realizable 

y(s) matrices and, if they can, to find G(s) whose 

minimal realizations are smallest and to assign poles 
wherever possible. It turns out to be convenient to 
answer these questions in terms of a reduced basis, whose 
definition is as follows. Let | 

, 1 


1. 

,17 

E10 
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-1.37 E10 s 2 

+4. 

.40 
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,06 
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E10 
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s 

+1. 

,18 

E10 
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s 

-1.07 

E10 

5 

s 

-1. 

,24 

E9 
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s 

+1.49 

E9 
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-1. 
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E9 
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s 

+1.99 

E9 
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s 

-2. 

,07 

E8 
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s 

+1.99 

E8 

8 

s 

-5. 

,55 

E6 

10 

s 

+6.93 

E6 

9 

s 

-4. 

,80 

E4 

11 

s 

+8.33 

E4 

10 

s 

-7. 

,09 

El 

12 

s 

+2.38 

E2 

11 

s 

+4. 

,5 

E-l 

13 

8 

-1.00 

E0 

12 

s 


Fig. 6. 

, Polynomial 

Matrix 


I 

In this section, we report on FORTRAN and PL/I 
Softwares developed to implement the MDP Algorithm on a 
digital computer, and our experience in the application 
of the software to the jet engine control problem de- 
scribed earlier in the paper. Both the programs use 
floating point arithmetic to implement the MDP Algorithm, 
considered over the field of real numbers. The FORTRAN 
version, using double precision arithmetic, affords 15 
digits of precision (decimal) on an IBM 370/158 com- 
puter. The PL/I version, using extended precision 
arithmetic, carries 33 significant digits. Our jet 
engine minimal design problem comes down to the question 
of determining the rank-seven kernel of a module morphism 
whose domain has rank nine, and whose matrix representa- 
tion In the usual basis contains thirteenth degree poly- 
nomials. In our experience, the principal difficulties 
arise fiom roundoff error occurring as a result'of finiLe 
representation of real numbers in the computer. 

There are two noteworthy features of the floating 
point KERPO (KERnel of a Polynomial Operator) software, 
pirst, it provides the user some control over the number 
bf digits considered significant during internal computer 
arithmetic. In actual problems, this appeared as the 
critical factor in obtaining acceptable solutions from 
the computer. Second, it performs a verification of the 
computed results up to four significant digits. Any 
discrepancy so pointed up, one attempts to rectify 
by varying the number of digits considered significant. ‘ 
In the case of the jet engine problem, after making 
several runs, we obtained an (apparently) acceptable 
solution from the PL/I version by setting the threshold 
for loss of significance near eleven digits. We can 
compare this solution with the known solution to the 
problem, represented by fig. 3. To di. this, we proceed 
as follows. 

The complete solution to the kernel problem appears 
in the form of seven elements in a rank nine module, 
which are the required reduced basis for the kernel. 
Represented in the usual manner, five of these contained 
polynomials of degree k^ one or less. It is Interesting 
to note that the existence of such elements can be pre- 
dicted by the following argument. We interpret fig. 3 


47.4 


Page 
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Appendix, on an IBM 370/158 computer. 

Starting with the numbers of table 1, together 
with the L and H matrices of table 2, we go through the 
calculations outlined In Section 3, and arrive at an 
exact rational-coefficient version of the kernel prob- 
lem of (22). By applying the MDP Algorithm, conceived 
now over the field Q of rational numbers, we are led 
finally to an exact reduced basis for the correspond- 
ing exact 2x9 matrix over Q[s). The seven basis 
elements turn out to contain polynomials of degree k^ 
equal to one and no polynomials of higher degree. Note 
that this means the floating point software missed at 
least two elements of first degree in the reduced basis. 
Rounded to fit in the available space, the seven basis 
elements obtained from exact software are indicated in 
table 4. , 

We would now like to compare the computer resources 
needed for the floating point calculation with those 
required for the exact calculation. 

In the floating point software, a sort of trial 
and error process was used to optimize the calculation 
by varying the threshold for loss of significance. 

Though this software did not reach a satisfactory an- 
swer for the jet engine problem, we have allotted from 
our experience about seven runs of two CPU minutes each 
to this calculation 


On the assumption that all solutions can be generated 
from the kernel basis, the logical conclusion is that 
'the two columns of (41) can be represented as a linear 
•combination of the five first degree elements in the re- 
duced kernel basis. Interestingly enough, the question 
of determining this transformation can itself be repre- 
sented as another kernel problem in polynomial modules. 
However, attempts to generate such a transformation 
turned out to be unsatisfactory. 

As an alternative approach to verifying the KERPO 
solution, we used two of the five first degree basis 
elements to realize a second order dynamical control 
scheme for the jet engine, along the lines of fig. 5. 
From fig. 5, we could then obtain a state description 
for the overall closed loon system, which we then com- 
pared with the corresponding optimal integral control 
scheme system of fig. 3. This comparison was based on 
the first few Markov parameters. Table 3 shows this com- 
parison for two of these parameters. 

Table 3 

Comparison of KERPO Results with 
Optimal Integral Control 


Each run occupied 400K bytes of 

t memory. 

Next consider the exact calculation. This soft- 
ware occupied 300K bytes of memory and- executed the 
jet engine calculation in 135 minutes CPU time. How- 
ever, the great majority of this time turns out to be 
consumed in Algorithm 2, which computes a reduced 
basis. This suggests strongly that more research on 

the reduction process a common one in the literature 

could have a considerably greater than average effect 

on practical applications of the method. Except for 
the reduction, the remaining part of the calculation is 
just about an order of magnitude away from being very 
reasonable; and improvements of that order can be ex- 
pected to occur in the near term, either through hard- 
ware or software advances. 

A comparison is made in table 5. Here it is seen 
Table 5 

Comparison of Floating Point and Exact Solutions 


We note that our solution appears to have identified the 
B and C matrices correctly, while it is in error so far 
as the A matrix is concerned. On Che basis of this 
evidence, we conjecture that roundoff error *ncurred in 
implementing the Euclidean division algorithm has the 
most serious impact on the correctness of the solution. 
This is because, intuitively, the effect of the A matrix 
In the state space corresponds to multiplication by 's' 
in the module. Since, in our case, the factors by which 
the matrix columns are multiplied are computed via the 
Bivision algorithm, we hypothesize this to be the source 
of the error. j 


I In order to solve the jet engine minimal design 
problem, then, one has the option of developing floating- 
point software which ha . increased sophistication or of 
^witching to softwares which permit exact rational cal- 
culations. The next section reports on the latter method 


that, on the average, the difference between floating 
point and exact softwares was about an order of magni- 
tude in computing time. 

For the exact solution, it is of interest also to 
examine integer sizes at various stages in the calcula- 
tion, Such a summary has been made in table 6. Note 

Table 6 

Integer Size During Exact Solution 


6. EXACT RATIONAL SOLUTION » 

t 

One way to avoid the difficulties of finite machine 
representation of real numbers is to consider the numbers 
pf table 1 as being rational instead of real numbers. 

It is then possible to get an exact solution to the jet 
engine problem, using softwares such as FORMAC or ALTRAN. 
These have the capability of rational and symbolic manip- 
ulation with an essentially unlimited degree of precision 
Naturally, as the calculation proceeds, one would expect 
the integer size to increase quite a bit. As a con- 
sequence, the storage requirements and computer time need- 
ed to manipulate these would also be substantial. In 
this section we give evidence as to the magnitude of 
these, especially to contrast with the requirements for 
the floating point calculations. This yields valuable 
insight into the tradeoffs involved in terms of computer 
usage needed to solve typical realistic jet engine con- 
trol problems from the polynomial approach. The results 
reported here are based upon FORMAC software written to 
Implement, in rational arithmetic, the procedure of the 


Stage of Computation 


j Resource 

/ Floating Point 

Exact 


f Algorithms 1 and 2 

Algorithm 1 

Algorithm 2 

Memory 

400K bytes 

300K bytes 

300K bytes 

CPU 

14 minutes (average) 

18 minutes 

117 minutes 


1. 

X 

State Matrices For Plant 

4 

2. 

Plant Transfer Function 

14 

3. 

Inverse of Closed Loop System 
<T) 

33 

4. 

Kernel Problem (2x9 matrix) 

45 

5. 

After Algorithm 1 

150 

6. 

20% through Algorithm 2 

270 

7. 

60% through Algorithm 2 

250 


Final Reduced Basis 

160 


Markov 

Parameter 

KERPO Solution 

Optimal Integral 
Control 

i CB 

-4.4968 

0.64815 

-4.4969 0.64814 

-0.13348 -0.17304 

-0.13348 -0.17305 

CAB 

2631.6 

1348 

70.451 3.544 

20693 

10890 

129.15 95.573 
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Table 4 

(Rounded) Reduced Basis from Exact Solution 


2.443E-3 

0.0 

0.0 

0.0 ; 

0.0 

0.0 

0.0 

1.601E-3 

0.0 

0.0 

0.0 

0.0 

0.0 

6.565E-4 

-3.824E-2 

1.125E-3s 

-7.199E-5 
-6. 649E-6s 

4.494E-4 

1.2B8E-53 

4.146E-4 

l,188E-5s 

4 . 146E-4 
5.151E-6s 

4.146E-4 

5.834E-6s 

0.0 

5.806E-4s 

-7.898E-3 
1.87 E-4s 

-6.41 E-4 
-1.331E-5s 

0.0 

3.653E-6S 

0.0 : 

3.322E-6s 

0.0 

-2.571E-6s 

0.0 

-J-768E-6S 

0.0 

8.467E-5s • 

j 

-1.526E-3 
6.54 E-4s 

9.685E-5 

-6.621E-5S 

2.893E-6 

-4.642E-5s 

0.0 1 
-4.423E-5s 

0.0 

-2.416E-5s 

0.0 

-1.847E-5s 

0.0 

-1.314E-4S 

0.0 

5.081E-4 

4.794E-4 

l,843E-4s 

4 .479 E-4 
1.713E-4s 

0.0 

-4.245E-5s 

0.0 

-3.076E-5s 

0.0 

1.279E-5S 

1 

0.0 

0.0 

0.0 

-2. 217E-6 

-3.491E-5 

+3.131E-6s 

0.0 

7.438E-6s 

0.0 

-1.421E-4s 

0.0 

0.0 

0.0 

0.0 

0.0 

1.622E-5 

0.0 

9.972E-5S 

2.935E-2 
8.69 E-4s 

2.949E-4 

-4.416E-4 

-4.068E-4 

! 

-2.25 E-4 

-2.434E-4 

0.0 

2.955E-4s 

that integer 

size before and 

after Algorithm 2 

is about 

in a minimal realization of N(s)D 

X (s) can be viewed as 


the same, while it nearly doubles during Algorithm 2. 

This also suggests that improvements in the efficiency 
Df Algorithm 2 may be possible. 

Finally, we summarize by commenting that ghe float- 
ing point software used on the order of 4 x 10° byte 
seconds of computing power, but eventually did not yield 
an acceptable solution. On the other hand, the exa^t 
rational software required on the order of 2.4 x 10 y 
byte seconds of computing resources and led to an exact 
solution. 

7. COMPENSATOR POLE ASSIGNMENT 

The exact rational software discussed in Section 6 
obtained the reduced basis, t^(s), 1 i 7, with 

'!<■>- K<.)j “ !> 

of table 4. From (20), where G^(s) is 2 x 2 and K^(s) 

Is 5 x 2, we see that the matrix N(s) must be 7 x 2 while 
D(8) is 2 x 2. Accordingly, 

r N(s). 

l D(i) J 

Ls a 9 x 2 matrix, which means from (32) that two kernel 
Elements * , i 


a linear combination of the determinants |d^(s)!cf (s) | 

Table. 4 makes it clear that |d| must have degree at 
least two; and so, since 

| : cf 6 1 - -1.4092974E-8s -4 . 7599915E-7 (47c 

> dj : dy | = -8.6656773E-8S 2 -2.9268875E-6s (47i 

|d, Id, | - 2.9060023E-8s. (47t 


-1.4092974E-8s -4.7599915E-7 
2 


jiust be chosen to effect a design. These elements (43) » 
will be linear combinations of the reduced basis elements 

U2>. If ^ 

ll x M : t 2 ( s )i 

has a linear dynamical interpretation as described in 
{he Appendix, then N(s)D”^(s) has a minimal realization 
whose state matrix has a characteristic polynomial 


Now let 


IdjOO I d 2 (s)|. 


£ ik d k (s >* f ik 


|d(s)| ■ | l f lk rf k (s>: 
k-1 1 K • 


c R, i 
7 


\d,:d-\ - 2.9060023E-8s, (47c) 

o * / 

with the polynomials in (47) serving as a basis for 
Rjfs], the R-subspace of Ms] consisting of polynomials 
of degree two or less, it is possible to construct an 
arbitrary polynomial j 

j |D(s)| - bjS 2 + b 2 s + b 3 , * (48) | 

jfor bj e R, i = 1,2,3 by forming an appropriate linear 
combination 

+ * 2 K: d 7l + 6 3 l rf 6i d 7'! (49) 

$ e R, i * 1,2,3. The B^'s are uniquely determined by 
tne b^'s. To complete a minimal design (48), it is 
only necessary to calculate f., and for k and 
j = 1,6,7. But certain results from tne exterior 
algebra, referenced in Section 8 , permit the calcula- 
tion of (f^, f 17 ) and « 21 , f 26 , f 27 ) as the basis 

of the kernel of the matrix [ 8 ^, - 82 * 8 ^]. Space pre- 
cludes a complete treatment of the theory, so we turn 
to the jet engine example. 

We make the selection 

|D(s)| - s 2 + 2s + 2, (50) : 

not so much because these dynamics are most desirable, 
but rather because the industry methods described in 
Section 3 could not be used to achieve (50) ir* a 
minimal ‘design. Thus, by solving this case, vc estab- 
lish potential superiority for MDP over existing in- 
dustry techniques. 

Starting then, with (50) and working backwards, we 


7 7 

I It 


j-1 J J 


,can calculate 8^, i * 
las well as (f^, f 26 . 


1,2,3 and thence (f ^ t *16’ *17 
f 0 _). These calculations were 


lk £ 2J l rf k (s) :‘V s) I * 


performed using exact arithmetic again. The results are 
presented here after rounding. First we obtain, 

8, - -4.202E6 (51a) 


k-1 j-1 * ■> 

by elementary properties of determinants. This shows 
that the characteristic polynomial of the state matrix 


8 2 - -1.154E7 
8, - -1.095E9 . 


(51a) 
(51b) 
(51c) ' 
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obtained from the basis for 
. A workable set of f. *s is 


the needed f 
ernel of [ii^» 

- MB, 


Outnut 2 


The f^.'s wepe, in our application, further scaled by a 
£actor J of 10 J to obtain the compensator pains as reason 
able numbers. This can be done without upsetting the 
compensator pole placement to be achieved. Eqs. (45) 
and (42) can then^be solved for t and t 0 of (43). The 
9x2 matrix (t *t 9 ), which represents our solution, is 


Fig. 7 Unit Step Responses. (a) Step 
on Inout 1 (b) Step on Input 2. 

As a final note in this section, it can be pointed 
out that the fixed poles in a compensator solution are 
the zeros of the greatest common divisor of the poly- 
nomials! Idj, I rfj I , i, j = 1,2,..., 7}. It is clear from 

the pairings (1,6), (1,7), (6,7) of our example that 
this GCD is 1, and thus that there are no fixed poles 
in the jet engine application. 


2.573 

1.687 

1.769s +1.183 
0.0202s -8.32 
•1.158s -1.607 
-3.076s 
0.744s 
1.622 

0.915s +6.582 


8. REMARKS 


8.1 Conclusions 

Considerable work has been done in the control 
systems area on polynomial design methods. Regardless 
which viewpoint one takes toward the definition of 
such problems, their solution is usually assumed to 
proceed according to algorithms of the type described 
jin the Appendix. Conceptually, this theory has achieved 
^considerable maturity, and so it seems appropriate to 
conduct an extensive case study of its application to a 
realistic problem. This is the reason for the jet 
engine control analyses carried out in t..is paper. 

The conclusions are generally positive in nature, 
though with some temporary limitations. On the positive 
iside, Sections 3 and 7 show that KDP is a problem rele- 
vant to the jet engine control industry and that the 
! MDP Algorithm offers a significant improvement in flex- 
ibility of design over existing algorithms in that in- 
dustry, The application problem detailed herein pro- 
vides a realistic and nontrivial test case for workers 
| in the area of computer solution of polynomial problems. 
A first limitation clearly occurs in Algorithm 2, which 
is a popular and well known theoretical algorithm. 

Both in terms of integer growth and relative CPU time, 
this reduction algorithm points to a need for further 
research. Following such an improvement, it would 
appear that the second limitation is overall CPU time 
for an exact solution. Though the cost of such time 
would be a small part of overall design cost, it appears 
desirable to reduce this tine by an order of magnitude. 
Since such a reduction seems to be a near-term possibil- 
ity by hardware or software advances, it would seem that 
polynomial methods may soon be ready to play a greater 
role in everyday practical design. 

8.2 Historical Remarks 


Matrix 


Elements 


0.1626 


The solution given in table 7 was verified by com- 
paring the corresponding Markov parameters for the 
closed loop systems of figs. 5 and 3. An exact calcula- 
tion comparing the first two Markov parameters, showed 
these to be identical for both systems. Another, non- 
exact calculation, which verified the first six Markov 
parameters, shoved agreement to four digits. The first 
two of these were listed in the second column of table 3 

J Step responses obtained from the closed loop system 
f fig. 5, using the numbers of table 7, are shown in 


The original stimulus for this work was the paper 
of Wang and Davison [1] in 1973, in which a minimal 
inverse system problem was solved. That work subse- 
quently led to the algorithm of Forney (2) phrased in 
rational vector spaces. Together, these works then led 
to the free-modular MDP Algorithm (3) which has been 
applied here. The jet engine application has been 
motivated by Michael and Farrar [4], whence arose our 
numerical data. A report on KERPO in double-precision 
FORTRAN has been presented (5), as has a more complete 
treatment of the pole assignment approach [6) in Section 
{7. Background reading on the algebraic aspects of the 
paper is available in [7]; and the exact proposition 
needed in Section 7 can be found in Chapter XV, Section 
B, Proposition 15 of (8). Further references to related 


A visual comparison of fig. 7,uitn similar plots 
Obtained for the optimal integral control system of 
fig. 3 showed them to be identical. Hence the latter 
$et is not included here. It might be interesting to 
examine the distribution of closed loop poles, which is 
l&ivcn below. 

-138.43 -4.47 + 0.986 i 

- 78.38 -1.678 + 0.238 i 

-0.136. 
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APPENDIX 


Let 


M : Fls]* 1 + F[s] P 


(A.l)^ 

be a morphism of free modules. In the appendix, we 
describe how a reduced basis for the kernel of M can be 
obtained and used to solve MDP. For a more complete 
discussion, the reader is referred to [3]. Solutions are 
pbtalncd in the form of the q x L matrix 


rncs y 

LD(s)_ 


(A. 2) 


U : FIs] 1 * - F[sT 


(A. 3), 


Step 1 . To the p x 
Identity matrix to form 


q matrix M, adjoin a q x q 


[fl 


<A.4)j 


Stop 2 . By elementary column operations, reduce 
(A. 4) lo the form 


E ll j 0 

E 21 • E 22 J 


(A. 5 


where has p rows, has no zero columns and is in an 
echelon torm. 

Step 3 . Then the columns of E.. are a basis for 
the image of M, and the columns of E 22 are a basis for 
the kernel of M (Ker M). 

Now, let b., i ■ 1,2,..., A be the columns of E^ 0 
obtained from Algorithm 1 as a basis for Ker M. Then 
by application of further unimcdular transformations, 
we can get an equivalent basis for Ker M which is re- 
duced in the sense of Section 4. Notice that we have 
introduced the notation b^ for elements of the basis 
before reduction, to avoid confusion with t-*, i « 1,2, 
-..A, which was assumed to be a reduced basis in Sec- 
tion 4. The algorithm below is used to reduce the ker- 
nel basis. However, the procedure is more general in 
nature and can be used to reduce A linearly independent 
elements in F[s]^ regardless of their origin. This one 
is typical of procedures described in the literature 
for doing these kinds of calculations. However, as has 
been pointed out in the paper, it is this part of the 
computation that consumes the major portion of computer 
time. Any research aimed at achieving efficiency in the 
reduction process is, therefore, the most likely to have 
a significant payoff in terras of making the MDP method 
of control system design tractable in the near term. 

Algorithm 2 

Write each b^, i = 1,2,. 

A 

b i - J b s j 

1 j *=0 1,J 

Where b. . c F^ and b. p ^0. We shall say the list 
^1* ^ 2 ***** is reduced if the matrix 


A in the manner 


(A. 6). 


Ib l.^ b 2 .t 2 

has rank A. Then, perform 


... b. . ] 
K * l K 


(A. 7) 


pver F[s], 

We change notation slightly by letting M be the 
fl x q matrix representing the morphism M. The technique 
s to choose q x q unimodular F[s]-mat rices to post- 
lultiply M. A matrix is unimodular if it has a non- . 
cro determinant that is an element of F. Mathemati- 
ally. 


Step 1 . If the list b^ , b 0 , • • • , of linearly 
independent elements is reduced, stop; otherwise, con- 
tinue. 

Step 2 . Determine field elements f. in F, l < i 
< A, which are not all zero and which satisfy 
A 

0. (A. 8) 


l Vi i 

i-1 1 lf4 i 


|ls unimodular if 1U| 0, e F. Such an operation is 

equivalent to a change of basis in F[s)^ and leads to a 
representation of M in the new basis. The following 
plementary column operations are examples of such trans- 
formations. The column operations are, (1) interchang- 
ing two columns of M; (2) adding an F[s ]-multiple of 
one column of M to another; (3) multiplication of a 
Column of M by a non-zero element of F. 

Given the p x q matrix M, the following algorithm 
leads to a basis for Ker M. The basis elements are re- 
presented in the usual manner. 

Algorithm 1 


Step 3 . For the set of integers i having f 
pon-zero, determine an i, denoted by i max , for which £, 

is a maximum, denoted by t , 

J max 

Step 4 . Perform the elementary column operation: 
Replace b^ by 
max 

! v i s< “ x ’* 1> 

i-i 1 1 

Return to Step 1. 

The question that remains is how the reduced basis 
may be used to obtain linear dynamical solutions to MDP. 
Let MDr take the form (A. 9) when stated over F[sJ. 


Now, in any solution 


"DM- 

lut lor 

[f] 


0. 


(A. 9) 
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£>f (A. 9), each of the £ columns will be contained in the 


kernel of M. All solutions pairs (N,U) can, thus, bo 
built up as linear combinations of elements in a basis 
for Ker M. Under what conditions will a solution pair 
(N,D) yield a minimal solution? 

Without loss of generality, we may assume that for 
lany candidate pair (N,D) the £ columns of 


b ± t 1 1 i 1 A 


(A. 15) 


for Ker M. j 

Step 2 . Apply Algoritlun 2 to the elements of (A. 15), 
to form a reduced basis. 


1 


[f] 


t £ , i < i <*. 


(A. 16) 


Step 3 . Express the reduced basis (A. 16) In this 
panner 


k i 

I k i * 
3-0 ° 


.3 


* ’ C i,k. 


iare 1 educed, because n they are not, a unimodular trans- 
formation V on F[s)' c , chosen according to Algorithm 2, 

jwJll produce an equivalent pair (N,D) such that the col- for i * 1,2,.,. A* Form the matrix 
limns of 


4 o 


m 


[are reduced and 

nd” 1 


NV(DV) 


-1 


kd" 1 . 


[t l,V x l 2,V. 2 


*.k. 


(A. 17 


(A. 18) 




(A. 10) 


£hen, ve make the following comments, offered without 
iproof . 

(1) K(s)D (s) can be realized by a linear dynamical 

system if ND“* is a matrix of proper rational func- 
tions. In such a situation, there exists a reali- 
zation A,B,C,E, all matrices over F, such that 


f the rank of the matrix formed from the last £ rows 
bf (A. 18) is net equal to £, stop; MDP has no solution; 
otherwise, continue. 

Step 4 . From the elements of (A. 16) in the reduced 
basis select £ elements 

t. i t. » • • • » t , 
with the properties 

(i) the rank of the matrix formed from the last £ 


G(s) - N(s)D _1 (s) 

- C(sI-A) _1 B + E. 


(A. 11) 


rows of 
lt ± k 


Equivalently, we also say that a pair (N,D) has a 
linear dynamical interpretation if the £ columns of 


W 


h u 1 


phese £ columns, when expressed as 

j 


pre such that the last £ rows of the matrix 

[^i f ^ • • • t * ) 

A,m^ 2, m 2 ^m^ 

have full rank. 


and 

■m 


are reduced and furthermore, letting the i^ 1 column be 


(ii) 


is equal to £; and 

£ 


A 


is a minimum. 


1 , 2 ,..., 


As a matter of fact, more solutions Co MDP may be 
j possible. Any elements t^, t ? ,...,t£ in Ker M, 

(A. 12) which admit a linear dynamical interpretation and achieve 

E he minimum order dynamics predicted in (ii) above, are 
solution to MDP throug \ the equations 


■i 

l hr 

3-0 i ’ J 


I s 0 


< A * 13 > and 


[t> 


G(s) 


Itj t 2 ... t^] 


N(s)D 1 (s). 


(A. 19) 


(A. 20) 


(A. 14) 


1 How, in the jet engine problem, A is 7 and £ is 2. The 


two needed columns of 


Being concerned with finding a realization with the 
Least order of dynamics, we state two more properties. 

(2) If the property in (1) is satisfied, then the deter- 
minant |D(s) | is related to the minimal realization, 
being an F-multiple of the corresponding charac- 
teristic polynomial |sI-A)|. Also, 

1(3) The columns of 




were generated in Section 7 to satisfy the pole place- 
ment requirement. 


[f] 


when expressed as in (A. 13), yield the number of 
dynamical elements in the minimal realization as 

£ 

r 

L V 

i-1 1 

With these notions, let t^, i - 1,2,..., A be the 
Reduced basis obtained from Algorithm 2. Then, MDP re- 
duces to generating t elements in Ker M which have a 
linear dynamical interpretation, with minimum order dyna- 
mics. For this, we can use the MDP Algorithm. 

I 

MDP Algorith m 

l Step 1 , Apply Algorithm 1 to obtain a basis 
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Abstract 

Present research efforts in the area of linear multivariable control systems in- 
clude activities which will probably reestablish frequency domain methods as fre- 
quently used tools for design. Two notable branches of this activity are polyno- 
mial methods and return-d if f crencc-deterninant methods. This paper sketches some 
features of these approaches, in the coutcxi of a numerical example from turbofan 
engine control. 


1. INTRODUCTION 

State variable methods for the design of linear 
multivariable control systems are well established 
as a major tool in the applications. Variants of 
the linear quadratic regulator theory are probably 
the most successful, with a variety of other tech- 
niques such as pole placement, decoupling, and 
geo: etric regulator theory also available. Even 
today, however, linear quadratic regulator theory 
still requires a somewhat indirect thought process, 
a feature it shares with many optimization methods; 
and much of the remaining technique is synthesis 
oriented instead of design oriented. 

Accordingly, some modern re-emergence of frequency 

domain thought has occurred especially for design 

Broadly depicted, this work involves polynomial 
methods and return-dif f ercnce-deterninant methods. 
This paper records certain studies of these ideas, 
on a conaian illustration from turbofan engine con- 
trol. Brevity precludes in-depth treatment; ve re- 
ly instead on the illustrations and the references. 

2. ILLUSTRATIVE PROBLEM 

The turbofan engine model chosen for t be Illustra- 
tions has two control input* — fuel flow and exhaust 
area* five states — fan turbine inlet tenperaturc, 
main burner pressure, fan speed, high compressor 
speed, and augmentor pressure, and two outputs — 
thrust and high turbine inlet temperature. In 
traditional (A,B,C,D) form, the state description 
( 1 ] is given by the matrices at the top of the 
following column, at a power lever angle of 47 *, 

For the sequel the corresponding matrix C(s), name- 
ly C(sl-A) is recorded. 

The design problem is to select compensators for 


57.095 3.613 -10.211 - 5-681 - 2.715 

19.832 -72.360 30.295 60.972 15.327 

.660 4.496 - 3,601 - .011 - 2.803 

1.326 2.313 - .809 - 3.032 - .821 

.832 .703 2.922 1.471 - 4.596 



C(s), in a loop under unity negative feedback of the 
plant outputs. Fast step responses with small over- 
shoot are of interest. 


3. POLYNOMIAL METHODS [2] 

Polynomial methods take advantage of the fact that 
action of the A-matrix and the s-variable are close- 
ly related in a nodule theoretic sense [3). Not yet 
well advanced computationally, polynomial methods 
nonetheless of ter considerable insight into system 
structure.* As is to he expected, they resemble the 
geometric methods in this regard. 
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C(r.) 


(.OlCs 5 ♦ .1453** - 92.05s 3 
~ yjC .% 2 + 29201s 4 95491) 


(-.0C6s 5 4 31.63s 4 4 3321.5s 3 


( .546s* 4 ?1.9s 4 2245s 3 
-1943. s 2 .166553 -12495) 

(-.033s 5 - .437s 1 * 4 63. 2s 3 


4 25500s 7 - 4 76068s 4 ?62?7) 4 l?03.3o 2 4 1742. 9s -3532.21 


b 5 4 140.7s 4 t 5337.6s 3 4 3c651s^ 4 119690s 4 133369 


— 1 the domain of the map represented 
by this matrix to determine seven 
"reduced basis" elements, shown be- 
low which serve to describe the 
kernel. From these, construction 
of K.(s) and K 0 (s) involves two 
linear combinations of these seven 
module elements, and standard reali- 
zation methodology. Using first, 
sixth, and seventh elements, and 
the assumptions 


V ■ V 


- 0, 


As an example, consider the selection of K (s) and 
K^(s) in Figure 1 in order to achieve a specified 
closed loop performance T(s). Such a rpecif ication 
is, of course, a nontrivial issue in its own right. 
A complete treatment of such a spec! f ication can 
be found in (2). Relying upon the algebraic inter- 
pretation of a transfer function as a pair of poly- 
nomials, such a design problem can be converted to 
a kernel calculation in R(s )-modules, where R(s] 
denotes polynomials in s with coefficients in the 
real number field K. Considerable manipulation 
must be carried out to sec up this kernel problem, 
which turns out to involve a 2x9 matrix of polyno- 
mials up to the thirteenth degree, as shown below. 


i 




] 


1 . 

17 

T 

E10 

3 

s 

- 1 . 

.37 

E10 

♦4. 

40 

E10 

4 

s 

-4, 

.39 

E10 

44. 

06 

E1C 

5 

s 

-4, 

.21 

E10 

41. 

18 

E10 

6 

s 

- 1 . 

.07 

E10 

- 1 . 

24 

L9 

7 

s 

41, 

.49 

E9 

- 1 . 

26 

E9 

8 

s 

41, 

.99 

E9 


07 

E 8 

9 

s 

41 

,99 

E 8 

-5. 

35 

E 6 

10 

s 

46, 

.93 

E 6 

-4. 

80 

E4 

U 

s 

48, 

.33 

E4 

-7. 

09 

El 

12 

s 

42, 

,38 

E2 

44. 

5 

E-l 

13 

s 

- 1 , 

,00 

E0 


10 

5 

11 

9 

12 


realizations can be found in the manner 
Compensator Realizations 


V 


ix 


Elements 



- 2.0 


-0.1626 



12.298 


0 



2.5732 


1.686 



-0.9369 


~ 65.039 



-3.2 


0.1 



1.092 


0 



-11.727 

-9.725 

2 . 8 S 8 

5.944 

0.824 

36.416 

3.28 

-13.66 - 

37.83 

9.147 

0.119 

0.777 

2.363 

9.973 

-3.806 

1.932 

0.022 

-1.265 

-3.36 

0.813 


Responses to unit steps in the two reference chan- 
nel arc shown in Figures 2 and 3. 

Solution of a problem by polynomial methods involves 
at this time nontrivial computational overhead, 
which is discussed in greater detail in [2]. It is 
likely, however, that advances in software and hard- 
ware will soon reduce this overhead. Advantages of 
the method include a finite enumeration of all solu- 
tions for a given T(s), and perhaps eventually a 
finite description of all possible performances. 


Sol utio n involv e s autonorphlc transf o rma t ions o n 4 . R F Ti^V-.nTrr^t>rvrr~nrT FRMI\ , A?n , S [4] 


T.443E-T - 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1.601E-3 

0.0 

0.0 

0.0 

0.0 

0.0 

6.565E-4 P rescnt computa- 

tional. situation for 

-3.624E-2 

-7.199E-5 

4.494E-4 

4.146E-4 

4 . 146F.-4 

4.146F.-4 

0-0 polynomial methods 

).3 25F.-3s 

-6.649E-6s 

1.283E-5s 

1.18SE-5s 

5.151E-6S 

5.834E-6S 

3.806E-4sj aa j ccs alternate fre- 

-7.698E-3 

-6.41 E-4 

0.0 

0.0 

0.0 

0.0 

q q quency domain approach- 

1.87 F.-4s 

-1 .331E-5s 

3.633-6s 

3.322E-6s 

-.2. 5 7. IE - 6 s 

-1.70SE-6s . 

8.467r-5s es of Interest. 

-1.526E-3 

9. 635E-5 

2.893F.-6 

0.0 

0.0 

0.0 

0.0 If ve set to zero 

6.54 E-4s 

-6.621E-5s 

-4.642E-5s 

-4.423E-5s 

-2.416E-5S 

-1.847E-5s 

-1.314E-4sa n( j denote by K, 

q q ve have the archetypal 

1 279E-5«J jn * t y negative feed- 
back preconpensat ion 

0.0 

5.0S1E-4 

4. 794K-4 
1.843F.-4s 

4.4794F.-4 

1.713E-4s 

0.0 

-4.245E-5s 

0.0 

-3.076E-5s 

0.0 

0.0 

0.0 

-2.21 7£-6 

-3.491E-5 

43.131E-6s 

0.0 

7.438E-6s 

0.0 problem. If K is 

-l*421E-4% ssucltc j to have state 

0.0 

0.0 

0.0 

0.0 

0.0 

1.622E-5 

0 0 description 

9. 972E-5^^ f D^) f then a con- 

7.935F.-2 
8.69 E-4s 

2.949E-4 

-4.416E-4 

-4.068E-4 

-2.25 E-4 

-2.434 E-4 

0.0 bined state descrip- 

2.955E-4^i on (A cr B c ,C c ,l) c > for 
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Fig. 2 

Closed Loop Response to Unit Step in 
Commanded Thrust; Polynomial Design 
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Fig. 1 

Compensation Scheme for Polynomial Design 
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Closed Loop Response to Unit Step in 
Commanded Temperature; Polynomial Desigi 


Nyquist Plots for Compensated System 
Direct Expansion 


Nyquist Plots for Compensated System 
Alpha Expansion 
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Fig. 6 

Nyquist Plot for dct(I*f€K) 
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Thrust 


Temperature 


Fig. / 

Same as Fig. 2 for Direct Approach Design 


CARDIAD Plot, Column 1, Uncompensated 


Seconds 


Thrust 


CARDIAD Plot, Column 2, Uncompensated 


Same as Fig. 3 for Direct Approach Design 

the Jcop can be obtained by an isomorphism on the 
product of the state spaces X and associated 
with the plant and compensator, respectively, pro- 
vided that the gain matrix DL> has no negative unit 
eigenvalues. For this situation, one has the 


important relationship that 

IhddJ |sI-a c |-|i+ck| | s i-a| IsI-aJ, 

upon which a Nyqulst study can be based. We refer 











expansion chosen for the retum-dlf fercncc-detor- 
minaut. The obvious expansion, shown for the pre- 
sent illustrative case, is 


and a less obvious, more recursive expansion in an 
WxN case is „ 


Design based upon Nyquist plots of |l+CK| is cade 
challenging by the intricate way in which the com- 
pensator K relates to the determinant. At present 
only introductory design studies based upon the 
expansions above have been made (4). An illustra- 
tion is the compensator 


r» 1000(s+l) f 

1 U s(s+200) I 

whicli was chosen by a cut-and-try method to increase 
the speed of response of the second output. Figures 
4 and 5 show the terms in the "obvious" and "a" ex- 
pansions for the compensated system, with Figure 6 
indicating the sun, exclusive of the unit terrs 5_n 
each expansion. Closed loop responses to reference 
steps in each channel are shown in Figures 7 and 8. 
Though the temperature response in Figure 8 is ac- 
ceptable, the thrust response in Figure 7 exhibits 
overshoot; and considerable interaction is evident. 


Fig. 11 

CARDIAD Plot, Column 1, < 


In current practice, plots such as Figure 6 tend to 
be the most useful. Design technique tends to focus 
upon rcducing^the interaction evident in these re- 
sults, which brings us to the next topic. 


5. CARD 1 AD— A DOMINANCE APPROACH (5] 


In making a Nyquist plot of the determinant of re- 
turn difference, H. H. Rosenbrock, [7] has estab- 
lished that |l+CK| encirclements can be counted as 
the algebraic sum of the encirclements of the diago- 
nal elements of return difference (1+CK) — provided 
that a condition of "dominance" holds on (I4-GK). 

Tills means, in our case, that the off-diagonal ele- 
ment in a column is smaller in magnitude than the 
diagonal clement, as a function of frequency (s=jw). 
Related to this stability oriented usage of the dom- 
inance Idea is a corresponding requirement on the 
loop transmission CK, which is used to help with de- 
coupling closed loop performance. 


Selection of K(s) for this latter purpose, so that 
G(ju)K(Jw) is dominant on its columns, has been 
widely studied for the case in which K(s) is re- 
stricted to he a constant matrix. Much less has 
been accomplished relative to the choice of a 
dynamic K(«). 


Fig. 12 

CARDIAD Plot, Column 2, Compensated System 

to such studies as rcturr.-di f f crencc-dc terminant 
methods, because of the presence of |l4CK| in a key 
role. 


A new technique for this purpose is the CARDTAD 


Construction of a Nyquist plot is related to the 



-.ot, acronyncd for Compensator Acceptability 
jgion for D lo gon a l dominance. Compensators 


K(ju) 



1 

OO+Jy^j*) 


* 2 (J<‘>)+Jy 2 (>>) 

l 


arc assumed, without loss of generality Tor pre- 
tompensation. A CAKDIAD plot for column one of 
the uncompensated system is shown in Figure 9. 

Each circle corresponds to a particular frequency 
ii> f and acceptable (x^.y^) pairs must be outside 
dashed circles at the frequency in question. Note 
that >' » 0 and x, suitably negative will be ac- 

ceptable for all frequencies. Figure 10 shows a 
CARIHAD plot for column two. Acceptable (x^,y 0 ) 
pairs must be inside solid circles at the frequen- 
cy in question. 


The simple compensator 


K(c) 


1 ,7s 4- .44 

.05s + 1 

-10 1 


achieves dominance at all frequencies in both col- 
umns, as can be seen in Figures 11 and 12, which 
consist only of solid circles each of which inclu- 
des the origin. 


Horc detailed information about an application of 
this method to design and simulation of a turbofan 
englt v 1 control can be found in [6]. 


6. DISCUSSION 


3. S. MacLnne and C. Birkhoff, Algebra . London: 
The Macmillan Co., 1967. 

4. P. W. Iloppner, "Hie Direct Approach to Compen- 
sation of Multivariable Jet Engine Models," 
Tech. Rcpt. EL-779, Notre Dame, Indiana, April 
25, 1977. 

5. R. M. Schafer, "A Graphical Approach to System 
Dominance," Tech. Rcpt. EE-778, Notre Dame, 
Indiana, March 26, 1977. 

6. * R. M. Schafer, R. R. Gejji, F. W. Hoppner, W. 

E. Longenbaker, and M. K. Sain, "Frequency 
Domain Compensation of a DYNGEN Turbo fan 
Engine Model," Proc. 1977 Joint Automatic Con- 
trol Conference, pp. 1013-1018. 

7. H. II. Rosenbrock, Compu ter-A i ded Control System 
Design . New York: Academic Press, 1974. 


Recent activities in frequency domain analysis and 
design of linear multivariable control systems sug- 
gest a certain resurgence of this viewpoint in use- 
ful new ways. Though somewhat limited by space 
constraints, we have tried to give a glimpse of 
some of these methods in the context of a numeri- 
cal model from the turbofan engine area. Focus 
has been on polynomial methods, which bear close 
resemblance to geometric control methods in an 
abstract algebraic sense, and upon methods re- 
lated to the determinant of return difference. 

The CARDIAD plot, a new dynamical ap roach to 
dominance, has been illustrated. 
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ABSTRACT 


This work explores an alternative to existing methods which are 
commonly used to design controls for jet engines. Whereas most modern 
designs implement piecewise-linear quadratic regulators, this represents 
an attempt to obtain a global nonlinear optimal control for a two-spool 
turbofan jet engine. 

A necessary starting point, therefore, is to have a good nonlinear 
model on which to perform the control studies. Unfortunately, the only 
accurate existing models of jet engines are (1) linear analytical models 
valid only for small regions, or (2) massive nonlinear, non-analytical 
computer programs which attempt to match experimental data. What is 
needed for this study is something which lies between these two extremes, 
i.e., a nonlinear, analytical model. 

A fifth order nonlinear model was developed in this study which cor- 
rectly models most of the qualitative behavior of the jet engine, but which 
fails to achieve strong numerical agreement with DYNGEN, a reliable non- 
analytical simulator. Several linear models were derived, both from the 
nonlinear analytical model, and also from DYNGEN. A time optimal control 
problem was formulated, subject to various constraints. Dynamic Program- 
ming theory and the Successive Approximations technique were explored, and 
applied to the problem of interest, while several improvements in the 
numerical programming were introduced. Analytical and numerical results 
were obtained for several models, both constrained and unconstrained. 
Finally, these results were tested on the two principal simulators, TYNGEN 
and the analytical nonlinear model. 

The study successfully achieved time optimal feedback control la*.s 
for various models of the two-spool turbofan jet engine. Furthermore, 
valuable insight into the nature of the problem was obtained, and much 



useful computer software was developed. However, an optimal control law 
obtained from any model can only be as good as the model itself. For 
this reason, more work is needed to develop a better nonlinear analytical 


model of the two-spool turbofan jet engine. 
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CHAPTER I 


INTRODUCTION 

This work explores an alternative to existing methods which are 
commonly used to design controls for jet engines. Whereas most modern 
designs implement piecewise-linear quadratic regulators, this represents 
an attempt to obtain a global nonlinear optimal control for a two-spool 
turbofan jet engine. 

A necessary starting point, therefore, is to have a good nonlinear 
model on which to perform the control studies. Unfortunately, the only 
accurate existing models of jet engines are (1) linear analytical models 
valid only for small regions, or (2) massive nonlinear, non-analytical 
computer programs which attempt to match experimental data. What is 
needed for this study is something which lies between these two extremes, 
i.e., a nonlinear, analytical model. 

Finally, after a suitable model(s) of the F-100-like jet engine is 
obtained, a time-optimal control can be calculated. This control will 
be determined subject to various constraints. It will be derived using 
Dynamic Programming and the Successive Approximations technique. 



CHAPTER II 



TWO SPOOL TURBOFAN JET ENGINE MODELS 
2.1 Introduction 

In this chapter, a hierarchy of models for a two spool turbofan jet 
engine is discussed. The configuration for this engine has been speci- 
fied by NASA Lewis Research Center personnel. A preliminary version of 
this work is given in reference [1] and has also been reported in [2]. 

The models have been classified as follows: 

Model 0. The actual jet engine (hypothetical). 

Model 1. The DYNGEN [3] simulation program, coded with data 

presumed to have been taken from experimental measure- 
ments on Model 0. This model solves 16 nonlinear 

* 

differential equations and uses data maps and thermo- 
dynamic tables which cannot be expressed analytically. 

Model 2. This model involves the primary thrust of this chapter, 
and is a 5th order nonlinear analytical model. It 
includes the 5 state differential equations which gov- 
ern the dynamic behavior of the system, along with 20 
algebraic equations which express the relationship 
between various engine variables. 

In addition to these nonlinear models, several linear models have 
been developed. Their original purpose was to provide an indirect method 
to compare Models 1 and 2. Subsequently, they also became important in 
the determination of a time-optimal control for the jet engine, when 
comparisons showed marked differences between Models 1 and 2. 

Model 1L5. This is a normalized 5th order linear model which is 
obtained numerically from Model 1, using the experi- 
mental DYGABCD [4] program of L. Geyser. 
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Model 1L3. 


Model 1L2. 


Model 2L5. 


Model 2L3. 


Model 2L2 . 


This is a normalized 3rd order linear model obtained 
by means of an order reduction performed on Model 1L5. 
This is the corresponding 2nd order reduction of 
Model 1L3. 

This is a normalized 5th order linear model obtained 
by taking partial derivatives of the analytical Model 2. 
This is a normalized 3rd order linear model obtained 
by means of an order reduction performed on Model 2L5. 
This is the corresponding 2nd order reduction of 
Model 2L3. 


2.2 Model 2 Development 

There are several purposes for the devel6pment of Model 2. First, 
it enables one to readily see the basic nonlinear relationships between 
the engine variables. This allows one to gain insight into their static 
and dynamic behavior. Second, it is fundamental that an analytical 
model be available for the application of optimal control techniques. 
Finally, linear models obtained by partial differentiation of this model 
tend to have more structure (zero entries in the ABCD matrices) than 
those obtained numerically. This in turn gives the linear control 
designer more insight. 

Model 2 was intended to be an approximation of Model 1, based on 
the specified engine configuration. Theoretical relationships developed 
in references 15], [6], and [7] were employed a starting point and 
certain simplifications suggested in [8] were used. In various situa- 
tions, least squares and exact fits were made to theoretical forms, and 
if a theoretical form was unavailable, polynomial, linear, and exponen- 
tial forms were used, whatever seemed to best fit the situation. 



If 4fdi 


It 



In most cases, the variables used in Model 2 correspond to those of 
Model 1. A letter key provides consistency among the variable names in 
the following manner: 

P a pressure 

T a temperature 

U a specific energy 

V a volume 

W a flow 

Similarly, numbers in the variable names identify engine locations as 
per figure 2.1. Table 2.1 is a list of all variables used. 


High-pressure turbine^ 

Ductburner*^ \ 

Combustor^ \ AHPFXT /-Low-pressure turbine 



Before delving into the details of the model, certain decisions had 
to be made regarding the choice of state variables and the order of the 
system. There is no general agreement as to what the order of a jet 
engine system is. It is a physical, not mathematical, entity and thus, 
every mathematical model is an approximation to the reality. Naturally, 
the higher the order of the model, the more accurate the approximation 
should be. The order that was selected (5th) was a function of the 
accuracy required by the control study to follow. This contrasts with 
the DYNGEN 16th order model, but is not an excessively low choice, for 
even first order models could yield reasonable results. 

The most obvious states to choose are the rotor speeds, (1) N and 
(2) N . The other selections were (3) the burner pressure, P, , a vari- 

r h 

able which is strongly affected by changes in fuel, WFB; (4) the burner 
internal energy, , a variable which is related by a constant to the 
burner temperature; and (5) the afterburner pressure, , a variable 
which is strongly affected by changes in the nozzle area, Ag. 

Table 2.2 gives a listing of the inputs, states and outputs. In 
actual existing control systems, inputs (1) and (2) are used, along with 
movable guide vanes mounted throughout the compressor and fan stages. 
These vanes cause changes in the air flow in a manner similar to the 
bleeds used in the model. 

Tables 2.3, 2.4, and 2.5 respectively are listings of the constants 
used in the model, the design values corresponding to the specified 
engine configuration, and the nonlinear state equations. Note that the 
state equations are formulated in teims of intermediate variables which 
have a very real physical interpretation. 

Table 2.6 is a listing of these nonlinear relationships existing 
between the state variables and the intermediate variables. Some of 



TABLE 2.1 

SYMBOLS FOR VARIABLES 


Symbol 


Variable Description 


8 

CNC 

CNF 

FG 

N c 

n f 

P CMAX 

P FMAX 

P 21 


21 


T, 


50 

r 55 


U 4 

WAC 

WAF 

WA3 


W 


CMAX 


AW. 


CMAX 
WFB 
w 

FMAX 

WG4 


WG50 

WG55 

WG7 

ZC 

ZF 


nozzle area 

corrected compressor rotor speed 
corrected fan rotor speed 
thrust 

compressor rotor 'need 
fan rotor speed 

compressor pressure ratio at surge 

fan pressure ratio at surge 

fan exit (compressor inlet) pressure 

compressor exit pressure 

combustor exit pressure 

afterburner exit pressure , • , 

fan exit (compressor inlet) temperature 

compressor exit temperature 

combustor exit temperature 

high pressure turbine exit temperature 

low pressure turbine exit temperature 

afterburner exit temperature 

combustor internal energy 

compressor airflow rate 

fan airflow rate 

airflow rate into combustor 

maximum compressor airflow rate 

correction term for maximum compressor airflow rate 

fuel flow rate into combustor 

maximum fan airflow rate 

gaseous flow rate out of combustor 

gaseous flow rate out of high pressure turbine 

gaseous flow rate out of low pressure turbine 

gaseous flow rate out of afterburner 

compressure surge margin 

fan surge margin 
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TABLE 2.2 

INPUT, STATE, AND OUTPUT VARIABLES 



Variable Description 

Symbol 

U 1 

fuel flow 

WFB 

u 2 

nozzle area 

A 8 

X 1 

compressor rotor speed 

"c 

x 2 

fan rotor speed 

n f 

x 3 

burner exit pressure 

P 4 

x 4 

afterburner exit pressure 

P 7 

X 5 

high pressure turbine inlet energy 

U 4 

y i 

thrust 

FG 

y 2 

high pressure turbine inlet temperature 

T 4 


TABLE 2.3 
CONSTANTS 


Symbol 

Description 

Value 

J (AJ) 

mechanical equivalent of heat 

778.26 

G 

force of gravity 

32.174049 

R(RA) 

gas constant 

.0252 

Y* 

ratio of specific heats 

1.4 

P 2 

fan inlet pressure 

518.668 

I C (PMIHP) 

high pressure rotor polar moment of inertia 

3.8 

i f (pmilp) 

low pressure rotor polar moment of inertia 

4.5 

V 

COMB 

combustor volume 

1.65 

V 

AFBN 

afterburner volume 

49.77 

CVMNOZ 

nozzle thrust coefficient 

.9494 

N C DESIGN 

high pressure rc tor design speed 

10070 

(XNHPDS) 



N DESIGN 

r 

low pressure rotor design speed 

9651 

(XNLPDS) 



N 

combustor efficiency 

20.71175 

C PC 

compressor specific pressure 

,24 

C PF 

fan specific pressure 

.24 

S,B 

combustor specific volume 

.20279 

C PHT 

high pressure turbine specific pressure 

.22589 

c 

PLT 

low pressure turbine specific pressure 

.27938 

<|) (PCBI.C) 

percent of compressor exit air bled for cooling 

.16 

a (PCBLOU) 

percent of bleed air which leaks into fanduct 

.208 

B(PCBLHP) 

percent of bleed air put into high pressure turbine 

.726 

y(pcbllp) 

percent of bleed air put into low pressure turbine 

.066 
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CNF 

T 

21 

CNC 


50 


21 


W, 


FMAX 


State if 




dN 

„ 2 

(1) 

C 

dt 

= (— ) 

* W 


dN 

^ 2 

(2) 

F 

dt 

.30 J 

\ i_n; 

s j 


dP, 

* 

(3) 

4 

dt 

- v RT [t 4 w 

COMB 


dP 7 

Ry*T 

(4) 


= v [WGi 

AFBN 


dU. 

o 

3 

(5) 

4 

dt 

_ VB 4 rr 

V p ■ 

COMB 4 


TABLE 2.4 

EQUILIBRIUM VALUES 


Value 

Variable 

2.75 

PFMAX 

2.948255 

WAF 

11899.1 

w 

CMAX 

9873.95 

aw cmax 

23.9299 

P CMAX 

586.467 

WAC 

2.55142 

WA3 

1.02310 

WG50 

742.957 

WG4 

.98730 

WG55 

1467.47 

T55 

2892.04 

T 

l 7 

2103.47 

WG7 

25.3522 

FC- 

2.9960 

ZC 

203.123 

ZF 


TABLE 2.5 
STATE EQUATIONS 


State Equation 


: [C PC WAC(r 21 - V + C PHT WG50(T 4 " V 
-:C pp MAF(T 2 - T a ) f C pLI WG55(T 50 - Tjj)] 

i3 + WFB - WG4] 

+ - WFB - WA3] 

L’{WG4 - WFB - WA3} + y*{T„WA3 - T.WG4 + T. (l+rOWFB}] 
4 3 4 4 
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Eq. H 


Equation 


(1) CNF = 


N DESIGN 9651 

r 


(2) T 21 = T 2 + 214.2732 CNF - 48(A g - 2.948255) 


(3) CNC = 


N C DESIGN /T 21 /T 2 10070 /T^/518.668 


(4) T- = T 01 + 743.2722 CNC - 68(A_ - 2.948255) 
J zl o 


< 5 > T 4 = V C VB 


(6) T 5u = . 727 T 4 


(7) P 3 = 1.05944 P 4 


(8) P = -6.20568 + .0129774 T ?1 - .0185376 P 3 


(9) W_, AV = 261.01 CNF - 63.196 
FMAX 


(!0) P pikIAX = 3.516739 CNF - .23561 


(11) WAF = W FMAX +28.502 1 - e 


- 2 - 313268 < P FMAX- P 21 ) 


(12) W_ MAV = 137.54 - 457.987 CNC + 564.325 CNC 2 - 188.113 CNC 3 


(13) AW cmax = 6.492 - 4.9749 CNC 

(14) P„ w< „ = 26.43184 - 89.C'i84 CNC + 109.7243 CNC 2 - 36.5756 CNC 3 

CMAX 


(15) WAC = 


>^ 2 ^ 518.668 


M CMAX + ^“cMAX* 1 “ 3 


- 3662 ' P CMAX - 










these can be readily observed in a DYNGEN listing, while others are far 
more obscure. Equations (9) through (15) are approximations to the fan 
and compressor block data maps. 


2 . 3 Lineari z ation and Order Reductions of Model 2_ 

Model 2L5 was obtained through a very tedious and time-consuming 
hand-calculated linearization. The partial derivative of each state 
equation and nonlinear function was calculated, and then combined to- 
gether to form a linear model. This linear model was then normalized 
as follows. 

Let A be n x n and let B be n x m. Then each state derivative may 
be vrritten 


x. 

X 


n 


j=i 


a u x j 


m 


I b . . u . 

A ij j 


(2.3-1) 


Let the values of x at the design point be denoted x, and denote the 
normalized state variable as x : 


x 


i 



(2.3-2) 


Similarly for the controls: 



Combine (2.3-2) and (2.3-3) with (2.3-1), 


n 


m 


(x 


i X i > = l a ii x i x i + E 

11 i=l 3 3 3 i=l 


b. . u. u. 
ij J i 


(2.3-3) 


(2.3-4) 


and simplify, resulting in 
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‘i = l (aii !j_ ) x 4 + l ( b j 
j=l 2 x. J j=l 


(2.3-5) 


Thus, elements of the normalized matrices are obtained by 


3 ij = a u 


(2.3-6) 


b.. = b.. J 
ij iJ T* 


(2.3-7) 


The normalized linear Model 2L5 is given in Table 2.7. The eigenvalues 
of Model 2L5 are (1) (2) -7 . 2264 + 1. 3913j , (4) -73.554, (5) -153.27, and 
(3) -343.11. The numbered eigenvalues can be associated with the state 
of like number, as they bear a loose resemblance with the diagonal terms. 
Note that all eigenvalues are negative, and the model is clearly stable. 

It seems quite reasonable that lower order models would be almost 
as accurate, suggested by the clear difference in the magnitudes of the 
eigenvalues. They will also be much easier to use to perform Dynamic 
Programming studies, saving much storage space and c.p.u. time. As 
mentioned above, the eigenvalues show that the states which will be 
eliminated as the order is decreased, are P^(3), then 11^(5), then P^(4). 

The method used to perform the order reduction is to first rearrange 
the states into partitions of "states to keep", XI, and "states to 


eliminate", X2. 


A 11 A 12 1 ill + M — 


A 21 A 22 J l 2 


(2.3-8) 


Now set the derivative of X^ equal to zero, since their dynamic behavior 


is to be eliminated. Thus, 


/ 
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A 21 X 1 + A 22 X 2 


B 2 u = 0 


(2.3-9) 


Solving for YL^, 

h - - A 22’ ltA 21 X 1 - B 2 < 2 ' 3 - 10 > 


is now replaced in equation (2.3-8) by its expression in (2.3-10), 
yielding 


X 1 ^ A ll ” A 12 A 22 A 21^ 


X 1 + (B 1 “ A 12 A 22 1 b 2 J - ( 2 - 3-11 ) 


Application of this method yields Models 2L3 and 2L2 per Tables 2.8 
and 2.9 respectively. Note that the most important eigenvalues have not 
changed significantly in the model reductions. 


2 . 4 Linearizations and Orde r Reductions of Model 1^ 

As previously mentioned, linear models obtained from Model 1 will 
be useful in comparing with those obtained from Model 2. In addition, it 
should yield a good model of the DYNGEN simulator in the area of the 
design point. 

The general method for obtaining numerical linearizations of Model 1 
is outlined in reference [10], including all the necessary program inputs. 
Additional insight into the selection of states for low order models is 
provided by DYGABCD. This stems from the identification technique used 
in DYGABCD, which is to perturb the inputs and states one at a time, and 
then measure the changes in each state derivative. A loose hierarchy of 
states in terms of their importance in the model is obtained by measur- 
ing how much each state perturbation affects the fan speed (which is 
certainly one of the most important states). A close resemblance with 
the choice of states for Model 2 occurs. 


Matrix Elements 


Eigenvalues 


-12.549 

.83305 

671.66 

-104.15 

50.953 


-1.6928 

-5.6135 

387.71 

21.135 

-55.855 


3.584 

1.645 

-392.67 

64.925 

-81.205 


.23786 

.14381 

-26.16 

-67.803 

-7.4745 


1.7011 

1.304 

160.53 

2.6314 

-105.74 


-343.11 

-153.27 

-73.554 

-7.2264 + 1.3913J 


0 

0 

1.2813 

0 

149.21 


1.4078 

.75817 

-122.31 

-48.928 

-3.092 


1.2138 

0 


TABLE 2.8 


1.461 

0 


MODEL 2L3 


Matrix Elements 

Eigenvalues 

-8.4226 

-1.2541 -.047955 

-72.576 

2.3957 

-5.9244 .0048531 

-7.1663 + 1 . 2441j 

-11.566 

56.677 -72.562 



3.406 .79722 
2.1241 .56161 
31.486 -64.49 


0 

0 1.461 


-.63298 

-.97907 -.01486 



0 1.2138 

1.072 .1598 
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TABLE 2.9 


MODEL 2L2 


Matrix 

Matrix 

Elements 

Eigenvalues 

A 

-8.415 

2.3949 

-1.2916 

-5.9206 

-7.1678 + 1. 2401j 

B 

3.3852 

2.1262 

.83984 

.5573 


C 

-.23288 

-.63061 

1.1412 

-.99068 


D 

.63396 

1.0656 

-.084644 

.17301 



TABLE 2.10 


MODEL 1L5 


Matrix 

Matrix Elements 

Eigenvalues 


-3.8 

-1.277 

2.067 

-1.152 

1.448 

-251.51 + 23. I47j 


2.748 

-5.39 

1.585 

-1.991 

1.071 

-96.366 

A 

377.9 

49.51 

-264.9 

86.807 

78.91 

-5.0491 + . 83858j 


31.26 

139.39 

-6.269 

-88.69 

27.83 



-176.5 

23.91 

-10.27 

-37.4 

-246.7 




-.00259 

.3553 






.2116 

-.31618 




B 


12. 54 

-13.774 






-.6201 

-99.3 






157.78 

6.84 




n 

-.8594 

-.1397 

.6672 

1.167 

-.1236 



.055591 

.00656 

-.001837 

.01354 

.85391 


D 


-.10277 

.90094 





-.013839 

.020856 
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TABLE 2.11 
MODEL 1L3 


Matrix Elements 


Matrix 


-2.4307 

3.8281 

2.4466 


-.70897 

-4.9579 

140.5 


-.81149 

-1.7235 

-94.982 


Eigenvalues 


-92.242 

-5.0644 + . 86196j 


1.395 

1.2585 

15.434 


.30875 

-.35303 

-98.208 


.034897 -.0083832 1.3734 

-.60014 .081353 -.12635 


-.023854 .86846 

.52351 .046147 


TABLE 2.12 
MODEL 1L2 


Matrix 


Matrix 

Elements 

-2.4516 

-1.91 

3.7857 

-7.5073 

1.2631 

1.1483 

.97844 

1.429 

.070274 

2.0232 

-.60339 

-.10555 

.19932 

-.55158 

.50298 

.17879 


Eigenvalues 


4.9795 + .914781 
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Table 2.10 details the results for Model 1L5 after normalization. 
Inspection shows that Models 1L5 and 2L5 do not closely match on an 
eleraent-by-eleraent basis, although there is not great disparity between 
the eigenvalues of the two models. 

Order reductions of Model 1L5 were performed by the same method 
detailed in section 2.3, rather than direct use of DYGABCD. It was 
felt that Model 1L5 was a reasonable approximation to the DYNGEN simu- 
lator at the design point, and it was desired not to rely heavily on an 
experimental program. Tables 2.11 and 2.12 list Models 1L3 and 1L2 
respectively. Again, the eigenvalues do not change appreciably after 
the order reduction. 

2 . 5 Comparison of Model Responses 

Concurrent with the development of Model 2 was the development of a 
computer program describing Model 2 (see Appendix). It employs an 
Euler integration and was used with a time step (DT) of .001, very 
suitable in light of the values of the eigenvalues of the linear models. 
DYNGEN was run with time steps of .01 and higher, for it employs a 
modified Euler technique [3] which allows larger increments to be used. 

In addition, the linear models were tested on program ABCD (see Appendix), 
with a Runga-Kutta integration. Figures 2.2, 2.3, and 2.4 show 

various time responses of Models 1 and 2. Unfortunately, the responses 
shown are the closest Models 1 and 2 came towards agreement. The linear 
responses are evidence that the linearization and order reductions were 
correctly calculated. 

Comparisons of the various linear models are given in figures 2.5, 
2.6, 2.7, and 2.8. There seems to be better agreement between models 
for the frequency responses involving the* high pressure turbine inlet 
temperature (T^). 
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State Space Trajectories - Model 1 vs. Model 2 
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Figure 2.9 is yet another comparison of Models 1 and 2. It shows 
the steady state equilibriums of both models as fuel is varied. It is 
not surprising that a change in fuel in Model 2 produces a corresponding 
change in steady state which is greater than Model 1 would produce. 

This follows since Model 2L2 is known to have higher eigenvalues than 
Model 1L2. Also shown are transients for step inputs between a low 
thrust point and a high thrust point (the design point) . 
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THE TIME OPTIMAL CONTROL PROBLEM 

3.1 Introduction 

In light of the disagreements between models found in Chapter II, 
no single model will he relied upon to determine a time optimal control 
law for the jet engine. The control problem Is to determine a controller 
which will drive a model from a low thrust equilibrium to the design 
point, in minimum time, and subject to certain constraints, as yet 
undetermined. In addition, it is desired that the controller be deter- 
mined by feedback control law, for the usual reasons: reduced sensitivity 

to plant variations; control over system stability; and regarding pro- 
gramming aspects, the ease at which a global solution can be obtained. 

3.2 Statement of the Problem 

The necessary first step is to reformulate the models as discrete 
time systems. Let 

x(t + At) = x(t) + At«f(x(t), u(t)) (3.2-1) 

represent the system with starting time k and terminal time N. It is 
understood that 

f(x(t), u(t)) = Ax(t) + B u(t) . (3.2-2) 

for linear models. Let x(k) be the starting state and let the terminal 
time N be defined as the first instant at which the system state reaches 
the designated target set S. All x(t) are eX, the state set. The per- 
formance index 

N-At 

J(x,u) « £ At (3.2-3) 

t-k 

t - k, k 4* At, » . . . , N— At 


S3 
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is to be minimized with u(t) s U, the control set, and u defined as the 
control sequence. 

ij = u(k), u(k + At) , . . .u(N-At) (3.2-4) 

Furthermore, the minimization is subject to hard constraints of the form 
g, (x(t), u(t» < c. (3.2-5) 

X X 

3.3 Constraint Determination 

The final step in a complete formulation of the control problem 
lies in the determination of the and of equation (3.2-5). There is 
a strong intuitive need for such constraints, for the physical jet 
engine has very real performance limitations. Brennan and Leake [8] 
have chosen turbine inlet temperature and surge margin as constraint 
variables in their studies of the drone engine, and similar constraints 
have been chosen for this study: (1) high pressure turbine inlet 

temperature (T4) ; (2) compressor surge margin (Z^) ; and (3) fan surge 
margin (Z ). The surge margin of a compressor or fan is defined as 

Z ^ ( P out/ ? in) - 1 (3.3-1) 

P - 1 

max 

If either the surge margin or the turbine temperature is too high, the 
constraints will be violated. By definition, let 

T 4 = g l “ft)) (3.3-2) 

Z = g„ (x(t), u(t) ) (3.3-3) 

c / 

Z F = s 3 (X(t)., u(t)) (3.3-4) 

The next step is to determine for each model. 

Model 2 presents no difficulty whatsoever since all three constraint 

variables are defined in the Chapter II development. It will be an easy 

matter to incorporate these equations into subsequent control tests. 
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The constraints are harder to determine for the linear models, and 
a starting point is needed. Control studies hy Basso and Leake [12] 
have used constraints which were strictly functions of the states. How- 
ever, such is not the case here. Simulations of Both Models 1 and 2 

show T, and Z to have very little steady state change over a wide range 
4 c 

of state space, yet step inputs elicit strong overshoots from both 
variables. Clearly the constraints must be functions of both the states 
and the inputs. 

Once again, DYGABCD was used to obtain linear expressions for the 
constraint variables. An order reduction was performed (per Chapter II) 
yielding the g^ for Model 1L2: 

T 4 - - .61059 x.,^ - .10759 x 2 + .50292 + .17689 u 2 (3.3-5) 

Z = - .20154 x. - .45813 x„ + .20423 u. + .14724 u_ (3.3-6) 

c 1 c i £ 

Z F = - .58229 x x + .46872 x 2 + .18877 u ± - .92545 u £ (3.3-7) 

Constraint functions were not determined for the other linear models, 
since Dynamic Programming solutions (see Chapter V) subject to constraints 
were only obtained using Model 1L2 and Model 2. 

The final task remaining is to determine reasonable values for c,. 
of the constraint equations. These will play a fundamental role in 
the optimal control solutions of Chapter V, for they are hard con- 
straints which will often affect the control chosen. After studying 
results of DYNGEN simulations, it was decided to use the following 
values: 

^ - .150 (3.3-8) 

c 2 « .105 (3.3-9) 

C 3 = * 


080 


(3.3-10) 



CHAPTER XV 


THE DYNAMIC PROGRAMMING METHOD 
4.1 Introduction 

It Has been pointed out in Chapter III that a feedback control law 
is desired, rather than an open loop control. Furthermore, the Dynamic 
Programming method has been extensively used in such situations to obtain 
numerical solutions. One of the more recent examples of its application 
is found in reference [12], where Basso and Leake have successfully ob- 
tained a feedback control law for a single spool turbojet engine. Use 
of Dynamic Programming methods to solve time optimal problems was shown 
to involve a successive approximations technique. 


4.2 Dynamic Programming Theory 

The basic applications of the Dynamic Programming method are fixed 

* 

time, free right end problems. Let 

x(t + At) = x(t) + f(x(t), u(t)) (4.2-1) 

with u(t) e U 


The starting time k is known, and the terminal time N is known. The 
target set is any x(N) e X. The object is to find 

(x) ** min (x,u) (4.2-2) 


for a given initial state x, where 

N-At 

X (x,u) - K(x(N)) + l L(x(t),u(t) t) 

t=k 


Rewriting: 


N-At 


J (x,u) = K(x(N) ) +L(x,u(k),k) + l L(x(t),u(t),t) 


t=k + At 


J k (x,u) = L(x,u(k),k) + J fc + ^ fc (x(k + At) ,u' 


(4.2-3) 


(4.2-4) 

(4.2-5) 
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The Principle of Optimality states that in order for the entire 
state trajectory to be optimal from k to N, it has to be optimal from 
k + At to N. Thus, equation (4.2-4) can be reformulated as 

J fc (x,u) = L(x,u(k),k) + V fc + At (x(k + At)) (4.2-6) 

which leads to 

V fc (x) = min {L(x,u(k),k) + + Afc (x(k + At))} (4.2-7) 

u(k) 

Since the minimization really only concerns u (k) , (u (k + At) .... u(N - At) 
are previously determined), u(k) can be defined as u. and Bellman's 
Equation results: 

V fc (x) = rain {L(x,u,k) + + Afc (x + f(x,u))} (4.2-S) 

u 

The boundary condition is 

V N (x) = K(x(N)) (4.2-9) 

These equations are necessary and sufficient for optimality. 

4.3 Successive Approximations Technique 

The task is to fit the time optimal problem (i.e., free time, fixed 
right end) , into a form which can utilize the basic Dynamic Programming 
method. This was developed by Leake, Liu, and Richardson in references 
[13] and [14], and later applied by Basso and Leake in [12]. 

As per [12], let V^ n (x) be any function such that V^ n (x) >_ V^x) 
and let v n (x,k) be a control law which results when performing the 
minimization. 

min [L(x,u,k) + V, n , .. (x + f(x,u))] (x,k) 4 S (4.3-1) 

u e U k t at 
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Tvfl 

It is shown in [14] that if V^‘ (x) is the performance index 

resulting from V n (x,k), then 

V fc (x) 5 .V k n+ 1 (x) <.V fc n (x) (4.3-2) 

and further that V^ n (x) converges monotonically to V^(x) in a finite 
number of steps, although each (x,k) may require a different number of 
steps. Thus, it is concluded in [12] tjiat 


V hl W = min [L(x,u,k) + V^(x + f(x,u))] 

* ueU - c*,k> 4 s 


(4.3-3) 


which very closely parallels Bellman's Equation. The only difference is 
that in the solution of the fixed time, free right end problem, equation 
(4.2-8) is relating two performance indices for the same state, but 
separated by At in time; whereas, equation (4.3-3) is relating two per- 
formance indices for the same state, and the same time k, one being a 
better approximation than the other. 

It now appears that the time optimal free time free right end prob- 
lem can be successfully solved, using the existing Dynamic Programming 
method. Indeed this is true for all practical purposes; however, there 
is a slight discrepancy between the successive approximation theory and 
its application to Dynamic Programming. To be specific, it is a fallacy . 
to conclude that equation (4.3-3) guarantees that (4.3-2) be true.' By 
definition, is the cost which results when applying V (x,k), until 

the target set S is reached, which is not equation (4.3-3). For example, 
let 


V ' °(x) - V 
k max 


(4.3-4) 


(x,k) 4 . S 


and 
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V, ° (x) = 0. . (4.3-5) 

(x,k) e S 

Then, if x is sufficiently far away from S, it is quite possible that 
there exists no v(x,k) e U which will enable the equation 

V^°(x + f(x,u)) = 0 (4.3-6) 

to be true, i.e. , the control could not cause the system to reach the 
target set in a time of At. Since equation (4.3-4) is true for all 
(x,k) e S, then 

V^fx) = min [L(x,u,k)] 4- V^Cx) (4.3.7 

ueU 

and equation (4.3-2) is no longer valid. In practical situations, how- 
ever, the method used in [12] and also used in. Chapter V of this 

study, using equation (4.3-3), will still converge. 

* 

A further simplification can be made when the control problem is 
time- independent, which is the case in this study (see section 3.2). 
Equation (4.3-3) simplifies to 

V n +I(x) _ m ^ n 4. V n (x + f(x,u))] (4.3-8) 

ueU (x,k) e S 


4.4 Technique Refinements 

One way of assuring that equation (4.3-2) will always be true is to 
replace (4.3-8) with 


V n ^ (x) = min 


V n (x) 


4 


min 

ueU [L(x,u) + V n (x + f (x,u))] (4.3-9) 


Rewriting this in terms of the problem as described in Chapter III, 

V n (x) 

min [At + V n (x + f (x,u))] 
ueU 


V U+ 1 (x) = min 


(4.3-10) 


A successive approximation problem allows still another departure 
from the basic Dynamic Programming problem (fixed time, free right end). 
Let us examine how equation (4.2-8), describing a fixed time problem, 
would be implemented on a computer. V^_^(x) would be calculated for 
all x e X, and stored; V jj_2At ^ would be calculated and stored, and so 
forth. Therefore, each Iteration has a specific time associated with it 
However, in the successive approximations technique, either all approxi- 
mations are concerned with the same time, or the problem' is time inde- 

pendent- In Basso and Leake [12], V (x) was calculated for all x e X 

•V .. - • ' ' • . ■' V, ■ ' 2 ' V 

and stored in .an array. Then V (x) was calculated, and after that had 

: '' ■. •' 2 x ’ • . 

been . completed far all x e X, V (x) replaced V (x) in the array, and so 
forth. It would be more efficient to immediately change each (x) to 
the just- calculated V' .(x) ih a- state-by-state maimer. In. reality then, 
the approximations- for V changes much more fapidly, for one does not 
wait uns.il the completion of the sweep through state space before using 
information derived during that sweep . In this manner, hew information 
becomes . available at a faster rate, speeding up the convergence to V (x) . 

Furthermore, i£ one starts the state space sweep at the target, and 
slowly moves away from the target, convergence will occur still faster. 
By starting near the target, one is testing controls for states, which 
can probably reach the target in a time of the order of At. Since V°(x) 
for the target is equal to zero, while guesses for V°(x) at other states 
must be made safely higher than the unknown solution, it is a benefit to 
start at a point where the information is the best, letting the informa- 
tion propagate outward to other states. Figure 4-1 shows the logic for 
this state search. This logic requires the target to lie at the center 
of the state space. '. 















TABLE 4.1 


A COMPARISON OF COSTS FOR "FAST UPDATE” SUCCESSIVE 
APPROXIMATIONS vs. THE “SLOW UPDATE” METHOD 

(V(x) . is listed first for each state) 










As proof of the numerical superiority of this "fast update” method, 

6 

Table 4.1 compares Dynamic Programming results for V (x) , one obtained 

through regular "do loop" sweeps through state space, the other by the 

"fast update" method. Both started with the same V°(x). Note that the 

superscript on the cost function no longer refers to the approximation 

number,, but serves as a record of how many sweeps through state space 

have been made. Xn the fast update of Table 4.1, there will have been 
2 

6p approximations made, where p is the number of discrete points for 
each state in tills second order system. Of course, c.p.u. time is 
virtually identical for both programs . 

Since the number of actual approximations in the fast update method 
is equal to the number of the sweep through state space times the number 
of points in state space, the finer the quantization, the more benefit 
is derived through use of the fast, update method. An alternative expla- 
nation is that the old method (do loops) makes you wait even longer 
before obtaining new information, when you increase the quantization of 
the state space. 

4.5 General Program Structure 

One of the first considerations is c.p.u. time. This is a function 
of the number of points in control space and state space, as well as the 
time increment At which is used. (It should also be mentioned that this 
refers to c.p.u. time on an IBM 370/158 computer.).. In this study, indi- 
vidual programs were limited to 14 minutes , 59 seconds, to avoid addi- 
tional job control language complications which occur for higher times . 
Thus, single Dynamic Programming solutions may be the result of several 
program runs. Near the end of the time limit, each job stores cost 
information on disk to be used by a subsequent job as a starting point , 


Input quantization information and time step 

ii 

— 

t 
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FIGURE ,4*2. Flow Chart of Dynamic Programming 














V° (x) . In addition, control information is stored, so that the optimal 
feedback control law will be easily accessed by the simulators in 
Chapter VI. 

Dynamic Programming is generally best-written in a somewhat ad hoc 
fashion. The number of subscripts in an array is dependent on the 
order of the system, and interpolation schemes will differ according to 
the dimension of state space. However, there still remains a basic 
structure to the program. Figure 4.2 shows a flow chart, while the 
actual program is contained in-, the Appendix. Note the absence of a 
T, do loop 11 for searching the state space. Also, the target cost is set 
to zero and left there, never allowing interpolation errors due to 
quantization to occur. When controls are tested for possible violation 
of constraints, the values of the "present state" x(t) are used. How- 
ever, the "future state" x(t + At) is used when testing whether or not a 
particular control takes the state outside of the state set X. The 
interpolation scheme is a standard method as used in reference [12] for 
two dimensions, and is analogously extended for third order models 
(Model 2L3 in particular). 


CHAPTER V 


OPTIMAL CONTROL LAWS 
5.1 Basic Problem Considerations 

In order to compare control studies of the linear models of 
Chapter II with studies of the nonlinear Model 2, linear affine models 
must be formulated. If the linear description of the system is 

x = Ax + Bu (5.1-1) 

■At A 

and the equilibrium values at the target are designated as x and u, then 
the linear affine system is 

x = (x - x) = A(x - x) + B(u - u) (5,1-2) 

Similarly the constraint variables become 

(y - y) = C(x - x) + D(u - u) (5.1-3) 

* 

where y is the equilibrium of the constraint variable. Since all the 

A A rt 

linear models found in Chapter II were normalized, x = u = y - 1. 

The time increment At (henceforth known as DT) for the linear 
models was selected based on the eigenvalues of each system. In all 
cases, DT = .01 seemed to be an acceptable choice, and convergence of 
the approximations did occur with this value. 

Quantization of the control and state spaces must be considered 
next. In general, one would like as fine a quantization as possible, 
but practical limitations on the cpu time will dictate a compromise. It 
is desirable that the quantization of the state space be small enough 
such that the program does not rely too heavily on interpolation. How- 
ever, if the state quantization becomes too small, DT must also be 
decreased. In other words, the amount by which a state can change in a 
time step DT will also have a bearing on the state quantization. For 



example, at a point in state space near the target, it is possible that 
the true optimal cost V(x) could be less than M, if the state quantiza- 
tion is too fine. 

The presence of constraints is important, for o'no desires small 
enough quantization to ascertain when the constraints are affecting the 

choice of control. If quantization is coarse, it may be much harder to 

* 

recognize that a control is riding a constraint. 

A big factor in an optimal control solution is the definition of 
the control set U, hot only as regards the quantization, but also tbe 
maximum and minimum values. These, of course, are chosen to reflect a 
true physical situation, and as such, it is expected that they influence 
the resultant control law. In these studies, the controls were limited 
such that 

0.5 5 WEB £1.4 (5.1-4) 

0.7<_A8 £ 1.2 (5.1-5). 

Again i these are normalized values. The state set, X, does not affect 
the solution for the states of interest, as long as these states are 
sufficiently far from the boundaries of X. 

5.2 Model 2L3 Unconstrained 

The basic choice of V°(x) for all models was 


V (x) = min 


V = .70 
max 


C l (x l ~ 1)2 + C 2 <X 2 " 1)2 + C 3 


with the c^ chosen such that V° (x) _> V(x) . Whether or not this condi- 
tion was satisfied was easily recognised by the success or failure of 
the V n (x) to converge to a solution. 
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Model 2L3 - Optimal Control Law (l Control) 
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FIGURE 5.2-A. Model 2L3 - Optimal Control Law (2 Controls) 
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FIGURE 5.2-B. Model 2L3 - Cost (2 Controls) 
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The program was originally developed using only one control, OTB. 
Subsequently, the A8 control was added but constraints were not yet 
considered. Figures 5.1 and 5.2 present the results, showing a single 
cross section of the three-dimensional state space, as defined by the 
plane x^CP^) = 1. 

It is interesting that the control law WFB(x) remains basically 
unchanged with the addition of the second input. However, the benefit 
derived by its addition is clearly evident, since the two-input system 
reaches the target roughly 10% faster than the one-input system, for any 
given state. 

Another significant result is that the optimal control solution 
for Model 2L3 is virtually the same for a given and x^, regardless 
of the choice of x^- The largest difference occurs at (x^ - 1, x^ = 1, 
x^ = 1) where V(x) — 0.0161 seconds, as opposed to V(x) = 0 at the target 
(x^ = 1, Xj, = 1, x 3 - 1) . The farther that x^ and are from the 
target, the smaller the difference becomes. This suggests that a second 
order model would be satisfactory to use, and, indeed, it comes as no 
surprise considering the eigenvalue information garnered in Chapter XX* 

Ey (afterburner pressure) reacts so quickly that it only slightly alters 
the cost when x^ is far from the target. 

5.3 Optimal Control Theory 

Henceforth, it is assumed that second order models of the jet 
engine are entirely suitable for this study. Considering the low order 
of the system, it was decided that analytical optimal control theory 
might provide good insight into the ultimate feedback controller 
solution, at the same time providing a means to check tha accuracy of 
the numerical program. The analytical approach is examined here. 


Let us first restate the general problem in continuous time: 


Minimize 


J(x ,u ,t) = K(x(t )) + L(x(t),u(t) ,t)dt (5.3-1) 

[t o ,ti J 

o 


such that 


x = fCx^Ujt) (5.3-2) 


where = the starting time 


t ^ = the time at which the target set S is reached 


X = x(t ) 
o o 


u(t) e II, the control set 


(x(t.j),t,) s S 


u r . . . = the continuous control over the interval from t to t 

rt 0 ,t] 


Note that the constraints g^(x(t) ,u(t) ) are excluded from the problem in 
this analytical study. 


For the problem of Interest, equations (5.3-1) and (5.3-2) can be 
restated as: 


Minimize 

J(jc o» u rt o ,ti» t) = 


fc i 


dt 


t- 

o 


(5.3-3) 


Such that 


x = A(x - f) + B(u — ' t) 


IT = {u(t) : u < u(t) < u . } 
mill — ■ — m a-TT: 


max: 


s = { (x(t^) , t^) : x = ?} 


(5,3-4) 


There are two principal analytical approaches to solve this problem, 
and both employ the Hamiltonian H defined as : 


E(x,p,u,t) = < p,f(x,u,t) > + p o L(x,u,t) 


(5.3-5) 


where p is known as the adjoint variable, and p equals 1 in this case. 
The first approach [14] states that if an infinum J exists for 
equation (5.3-3), 


J (x ,t ) = inf J(x ,t ) 
o o o o 


VL 




(5.3-6) 


then, it solves the Hamilton- Jacobi equation 

8J*(x o ,t Q ) + H*(x, 3J* T (x o ,t Q ) , t) 


- 0 


3t 


ax 

(x o ,t ) * S 


(5.3.7) 


(5.3-8) 


i 

; 


r ' 
1 ■ 

I 

j 


I 



J*(x ,t ) =* K(x ,t ) for (x ,£ ) e S 

O Q O O O O 

A 

Solution of equation (5.3-7) will yield a control v (x, 3J (x o ,t Q ) , t). 

3x 

If this control does carry (x(t o ),t Q ) to S, the target set t then the 

* 

control is optimal. Often the control law reduces to v (x) , i.e., a 
pure state feedback control law. This is similar to Dynamic Programming, 
and, in fact. Bellman’s equation (4.2-8) is actually a discretization of 
the Hamilton- Jacobi equation (5.3-7). 

In this case, an easier analytical approach is through application 
of the Minimum Principle [ 93 and the Hamilton Canonical Equations. 
Pontryagin’s Minimum Principle states the Hamiltonian must be a minimum 
as a necessary condition for optimality, i.e., 

HCx*,p*,u*,t) £ H(x,p,u,t) , (5.3-9) 

¥■ u e U 

& 

where u is the optimal control. The Canonical Equations state that 

p = - 3H~ (5.3-10) 

3x 

x = 3H T (5.3-11) 

5p 


such that 


< p,dx > - Hdt — p dK 

fc l 


0 


(5.3-12) 


where the differentials of equation (5.3-12) are consistent with the 
problem constraints. Equation (5.3-12) is known as the Transversality 
Condition [14]. 

For the problem of Interest, 


H = < p,Ax + Bu > + 1 - < x > + < B^p, u > + 1 (5.3-13) 


where p,x, and u are vectors. To minimize (5.3-13), u must lie at ah 
extreme point bf its control set, depending on the' signs of (B p)^ and 
(B T p)„ i.e. , 


* 

u = 


(1) 

U- 

1 max 

u 2 

max 

if 

(b t p ) 1 - - 

(B T p) 2 =* 

(2) 

1 min 

u 2 

max . 

'if 

(b T p) 1 - + 

(b T p) 2 = 

(3) 

U- . 

1 min 

u 2 

min 

if 


(b T p) 2 = 

(4) 

1 max 

U 2 

min 

if 

(^>1 = - 

(^>2 = 


(5.3-14) 


Thus we obtain the usual "bang-bang” solution, so characteristic, of many 
time optimal control problems. It is now certain that the controls will 
ride the boundaries of the control set U. The problem lies in determin- 
ing what control is applied when, and for how long. This leads one to a 


In reference [11] , Fontryagin shows a method for such an analysis. 
Equation (5.3-14) shows that if the trajectory of p (the adjoint system) 
is known, then the switchings are known,. Thus Pontryagin delyes into an 
analysis of the adjoint system, employing various transformations and 
translations to obtain 


Pl = - Xp 1 - w? 2 

k “-WPl - Xp 2 


(5.3-15) 

(5.3-16) 


where X and p describe the complex eigenvalues X + pj as per Chapter II, 
Further analysis of equations (5.3-15) and (5.3-16) yields the time 
intervals for which each control will be at 


5.4 Model 2L2 Unconstrained 

For Model 2L2, unconstrained, controls (1) and (3) as per equation 
(5.3-14) have time intervals of approximately 1.27 seconds, while con- 
trols (2) and (4) have time intervals of 0.00149 seconds. From this 
information, state space trajectories can be constructed. Since the 
solutions of (5.3-15) and (5.3-16) are basically sinusoids, the sequence 
of controls will begin at a particular control, depending on its loca- 
tion in state space, and follow the control sequence in numerical order 
Cand repeating) for the specified time intervals, until the target is 
reached. This, naturally means that the first and last intervals may 
be shorter than the others. 

Let us first construct trajectories for a number ot initial points 
in state space, given that only a single control is applied. As shown 
in Figures 5.3 through 5.6, each set of trajectories has a point of 
singularity, which is the equilibrium point of the system with the given 
control applied. Each "x" mark shows a time interval of 0.1 second. As 
the figures show, there exists a single trajectory for each control 
which will pass through the target. Let us start at the target, and 
reconstruct the trajectory for each control, going backwards in time 
(Figure 5.7). Considering the optimal control knowledge embodied in 
equation (5.3-14), it is clear that the final stage of any optimal tra- 
jectory must necessarily follow part or all of one of these arcs, in 
order to reach the target. Furthermore, if the last stage of an optimal 
trajectory follows the trajectory due to control i to the target, then 
it necessarily was moving in accordance with control law i-1 prior to 
the switching. In this way, figure 5.8 can he constructed. There will 
occur at most two switchings in the optimal control law for the area of 
state space which is of interest in this problem. 
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FIGURE 5.9. Model 2L2 - Optimal Control Synthesis 






5h 


By plotting the four trajectory systems on one graph, a composite 
picture is obtained, showing the optimal trajectories for the entire 
state space. This is shown in Figure 5.9. Due to the relatively short 
time interval for which controls (2) and (4) are applied, they have a 
negligible effect on the solution. In fact, for the scales used in 

act 

figure 5.9, these control regions do hot appear. 

As an example, suppose the initial state is (1.05, 0.75). Control 
(1) would be applied for approximately 0.2 seconds, then control (2) 
for 0.00149 seconds, and finally control (3) for approximately 0.06 
seconds. Control (2) could have been eliminated for all practical pur- 
poses, In fact, that is precisely Ttfhat occurs in the Dynamic Programming 
for the problem. 

The Dynamic Programming results for Model 2L2, unconstrained, are 
shown in Figure 5.10. In general, there seems to he good agreement be- 
tween the analytical study and the numerical results when considering 
the amount of time necessary to reach the target, V(x). The control 
laws, however, are not exactly the same. While the Dynamic Programming 
results show controls (1) and (3) to be optimum in the same areas (for 
the most part) as the analytical results, the boundary areas between the 
control regions do not agree as well as had been desired 

This is accounted for by a simple explanation. For some states, 
there exist control laws (in the Dynamic Programming results) which are 
not theoretically optimal. But when these control laws are applied, 
they yield costs which are so close to the optimal cost that, in a 
numerical study subject to interpolation error, a non-optimal cost with 
a corresponding non-optimal control law may be chosen over an optimal 
cost with a corresponding optimal control law. Thus, although the 
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difference between controls (1) and (3) is large, (where one would be 
optimal and the other non-optimal) , application of each at the particu- 
lar state (in the vicinity of a boundary) will yield costs which are 
nearly equal. 

5.5 Model XL2 Unconstrained 

Switching point analysis (per Section 5.3); for Model 112 shows that 
controls (1) and (3) have maximum time intervals of 1.22 seconds, while 
controls (2) and (4) have maximum intervals of 0.49 seconds. In this 
case, a maximum of one switching may occur for the state space areas of 
‘interest. Figures 5.11 through 5.14 show the various trajectories for 
each control. In this case, the singularity points are not as close 
together as they were for Model 212, and thus the different controls 
cause much different trajectories to occur. Recall that in Section 5.4, 
the close proximity of two singularities revealed that two controls both 
had the same general effect on the state trajectory, but one control 
was always slightly better for a much larger area of the state space. 

The composite effect is shown in Figure 5.15. These results are 
dramatically different from those of Model 112, which is not surprising 
since the relative magnitudes of elements in the B matrices of the two 
models is very different. 

f 

Although the results are not included in this work, it is a 
relatively easy job to devise a controller based strictly on these 
analytical results. Simple knowledge of several points along each of 
the four main trajectories (which completely define the control regions) 

will allow one to fit a curve to each boundary. .-. In this, way, the 

• • ■' * 

control is applied according to the region of state space, continuously 
testing, as the trajectory moves throughout state space. This would be 
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the method utilized in a controller simulation, if constraints did not 
need to be considered. 

The Dynamic Programming results for Model 1L2, unconstrained, are 
shown in Figure 5.16. Again, while there is excellent agreement on V(x) 
between the analytical and the numerical results, the optimal control 
laws do not precisely agree, for the same reasons as mentioned in 
Section 5.4. 

5 . 6 Model 1L2 Constrained 

Analytical results are not possible when the constraints, as 
developed in Chapter III, are considered. Basically, these state- 
control constraints may be interpreted as control constraints which 
vary as a function of the state. Furthermore, it is not clear which 
constraints affect the control law at a given state merely by studying 
the trajectory of an optimal solution. If the constraints were functions 
of the state only (as was the case in reference [12]), the optimal tra- 
jectory would easily reveal when the control was riding a state constraint . 

The Dynamic Programming results for Model 1L2, constrained, are 
presented in Figure 5.17, and the effects of the constraints are seen in 
Figure 5.18. Mote that each of the constraints has an effect on the 
optimal feedback control law for some area of state space. In fact, 
there is only a small region of state space where the constraints do not 
affect the solution. The main impact of these constraints is that the 
control law no longer even resembles a bang-bang controller, but instead 
is a continuously changing function. The more finely quantized that the 
control set U becomes, the smoother the control law will be. 
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5 . 7 Model 2 Unconstrained 

After sufficiently analyzing the linear models, it now remains to 
study the nonlinear Model 2. No attempt will be made here at a non- 
linear system analysis, which would be more complicated than the linear 
analysis outlined in Section 5.3. Furthermore, while the Dynamic 
Px'ogramming theory remains unchanged, the actual programming gets much 
more complicated. 

If Model 2 could be embodied in 5 state equations, each a function 
of state variables and inputs, there still would be no significant 
difference in programming difficulty from the linear models. Unfor- 
tunately, Model 2 consists of many equations involving many intermediate 
variables, and the number of operations required for each point in state 
space is dramatically increased. In fact, a fifth order solution would 
require a prohibitive amount of c.p.u. time, and is automatically ruled 
out. A method must now be found to reduce the order of the system, 
preferably to second order. Ideally, one desires to set the derivatives 
of the unwanted state variables to zero, just as was done in Sections 2,3 
and 2.4 with the linear models. However, it is impossible to obtain a 
closed form solution for all the intermediate variables, with state 
variables 3,4, and 5 eliminated. This leads us to consider iterative 
solutions. 

If equations 1 through 25 of Table 2.6 are to be solved (which is 
necessary to evaluate the state derivatives), then values for P^, and 
P_, must somehow be determined. Recall that Dynamic Programming involves 
a determination of x(t + At) for a given x and a given u. The first 
attempt, then, was to supply an initial guess for the eliminated states 
(P,, U . , and P 7 ) , holding N , N , and the controls fixed, and iterate 



until a steady-state solution was reached. 

This method failed for several reasons. If one undergoes this 
iteration process for each state and each control on every successive 
approximation, c.p.u. time is extremely high. Alternatively, if one 
stores all the steady state values for P_, , P^, and for each state 
and each control after the first approximation (eliminating the need to 
iterate on subsequent approximations) a prohibitively high amount of 
memory is required. Furthermore, if an initial guess for the eliminated 
states is not close to the actual steady-state solution, instabilities 
will occur and the system blows up. All of which requires us to look 
for another solution. 

As a compromise to the problems encountered in determining values 
for the eliminated states, linear approximations are obtained from 
Model 2L5 and an order reduction is performed. This eliminates any 
instability problems and also drastically reduces c.p.u. time. The 
resulting equations are 

P 4 = 1.4663 +.53032 x £ +.40998 ^ -.18155 u 2 (5.7-1) 

P ? - -.15936 +.78107 x 2 +.43393 -.88875 u 2 (5.7-2) 

U 4 = -.63063 x 1 -.99071 x 2 +1.0656 ^ +.17300 u 2 (5. 7-3) 

These equations are then converted to linear affine for utilization by 

the program. 

The Dynamic Programming computer progi*am for Model 2 is found in the 
Appendix. It is divided into four subroutines: (1) the main program, 

which is basically the Dynamic Programming method as outlined in Figure 
4.2, along with the constants used in Model 2 per Table .2,3; (2) the 
static relations of Model 2 per Table 2.6; (3) the dynamic relations of 
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FIGURE 5.19-A. Model 2 (Unconstrained) - Optimal Control Law 


! 


-3 



Vi 

1 1 




0.90 

0,92 

0*94 

0,96 

0,98 

:i 
> ! 




* 

* ' 

* 

* 

* 



1.10 

* 

0.2420 

0,2332 

0.2246 

0.2160 

0,2074 

*. i 


1.06 

* 

0.2245 

0.2146 

0.2052 . 

0.1956 

' 0,1859 

£ | 


1,06 

* 

0,2054 

0,1946 

0.1835 

0.1730 

‘ 0,1613 

i 7 


1.04 

♦ 

0.1344 

0.1723 . 

r 0 .1602 

0.1471 

0*1322 

i . 

r 1 


1,02 


0.1611 

0,1475 

0.1531 

0.1162 

0.0954 

3- 
J ' 


Nj , 1,00 

* 

0,1350 

0.1196 

0.1006 

0.0769 

. 0.0453 

t ; 


0.99 

* 

0.1069 

0,0923 

0.0734 . 

0 . 0 b 33 

0.0391 

r ! 

• 

0.96 

* 

0,0935 

0,0007 

0 1 0703 * 

0.0679 

0.0812 


* 0.0951 0.0912 0.0934 0,1037 . 0,1146 

0,92 .* 0.1131 0,1166 0,1246 0,1327 0,1403 

0.90 * 0.1363 0,1440 0,1500 0.1556 , 0,1609 

* ♦ # # $ 
0,90 0.92 0.94 0,96 0,98 

FIGtJRE % 5.19-B. Model 2 ( Uncons trained) - Cost 



1,00 

1,02 

1,04 

1*06 

1.08 

1*10 




* 

* 

* 

* 

* 



0,1986 

0,1095 

0,1795 

0.1685 

0,1559 

. 0,1416 

* 

1,10 

0,1755 

0.1642 

0,1515 

0*1366 

0,1193 • 

. 0.1056 

* 

1,08 

0.1483 

0,1334 

0.1156 

0,0945 

0,0832 1 

* 0*0002 

* 

1,06 

0,1144 

- 0 * U 929 

0,0670 - 

0, 0620 

0.0668 

' 0.0760 


■ 1,04 

0,0689 

0 . S 363 

0.0444 

0 . 0601 - 

0.0772 

0,0920 

*- 

1.02 

0.0 

0 . 03*3 

0.0629 

0,0833 

0,0997 

' 0,1120 

* 

1.00 

0,0572 

0,0752 

0,0950 

0 .1009 

0.1205 

0,1003 

* 

0,98 

0,0960 

> 0,1092 

0,1205 

0.1303 

0,1389 

0,1471 

# 

0,96 

0,1247 

' 0,1336 . 

0.1415 

0,1487 

0,1553 

0.1625 

* 

U .94 

0.1473 

0.1536 

0.1594 

0,1640 

0.1702 

0,1767 

* 

0,92 

0,1658 

0,1704 : 

0.1748 

0,1792 

0.1839 

0.1906 


0.90 

* 

* 

* 


* 

t 



1,00 

1*02 

- 1.04 

1*06 

1.08 

1*10 




-3 

ro 


73 


Model 2 per Table 2.5, minus the eliminated state variables; and (4) 
equations for determination of P^, P^ and U^, per equations (5,7-1), 
(5.7-2) and (5.7-3). Note that subroutines XSLOW and M0DEL2 utilize the 
unnormalized system, while the main program and XFAST utilize the 
normalized system. Thus conversions from one system to the other are 
made at several points in the program. 

It is also important to note that, even though c.p.u. time has been 
cut as much as possible, it still takes up to five times longer to obtain 
a Dynamic Programming Solution for Model 2 than for the linear models. 

For this reason, it is important to use as much c.p.u. time as possible 
per job run, but still leaving enough time to insure that the results are 
stored on disk before the allotted program time limit is exceeded. The 
program itself insures that the results are safely stored, by measuring 
how much c.p.u. time is required for the first successive approximation, 
and then using that information to decide when to write the results on 
disk. 

The first results presented in Figure 5.19 are for Model 2, uncon- 
strained, and are normalized values. The control law is similar to the 
solution for Model 2L2, unconstrained, but the cost is less than the 
Model 2L2 cost for most points in state space. It is somewhat surprising 
that nozzle area does not ride the limits of U at several states, 

5. 8 Model 2 Constrained 

Figure 5.20 shows the Dynamic Programming results for Model 2 with 
the constraint limits as specified in equations (3.3-8) , (3.3-9) , and 
(3.3-10), and using equations 23, 24, and 25 of Table 2.6. While the 
control can actually only ride one constraint at a time, there is a 
large area of state space where both T^ and 2 ^ are very near their 
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FIGURE 5,20-A. Model 2 (Constrained) - Optimal Control Lav 
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respective limits, and thus they both are shown as affecting the 

solution in Figure 5,21. If only the constraint closest to its limit 

is shown, it would result in at least 10 smaller regions. In fact, it 

is unrealistic to show these smaller regions, since the somewhat coarse 

quantization of U used in this study will have a strong effect on which 

constraint the controls are rising. The result is an almost random 

choice as to which constraint (Z or T,) that the controls ride (in the 

c 4 

T,-Z region). Note also that the constraints affect the solution much 

*T C 

differently than they did for Model 1L2 in Figure 5-18. 

This Dynamic Programming solution requires approximately 90 minutes 
of c-p.u. fcime s using a 225-point state space search 3 and a 209-point 

r 

control space search* 


CHAPTER VI 


CONTROLLER SIMULATIONS 

6.1 Introduction 

Now that feedback control laws have been obtained for several 
models, both constrained and unconstrained, any state (N^, N^) can be 
driven to the specified target in (approximately) minimum time. Thus, 
initial starting states can be chosen which reflect low-thrust condi- 
tions for both Model 1 and Model 2. The simulators read the feedback 
control law from its disk storage and change the fuel flow and the 
nozzle area accordingly. This will necessarily involve an interpolation 
scheme to decide where the trajectory lies in state space, and is accom- 

o 

plished using the same scheme as employed in the Dynamic Programming. 

It is desirable that the quantisation of the control law be as fine as 
possible for these purposes, but limits on c.p.u. time once again lead 
to a compromise. 

Interpolation will often lead to error when you are interpolating 
in a region of state space where the control laws change abruptly. 
Obviously, the optimal control which is desired is either one extreme 
or the other, and not something in-between. In this case, the interpola- 
tion scheme could be overridden by an analytical test, similar to that 
which x?as mentioned in Section 5.5. Such considerations have not been 
implemented in this study. 

The addition of constraints, and their resultant "smoothing" effect 
on the original "bang-bang" control law, will reduce the number of abrupt 
changes in the control law. Hence, the interpolation scheme is expected 
to be more reliable when applying a control law which was derived sub- 
ject to constraints. 
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There is another very fundamental consideration for the controller 
simulations. While the optimal control law' at some arbitrary state 
(x^, x^) is quite possibly the same, or nearly the same as that for 
state (x^ + e, + <5), (s,S small) , this is not at all true near the 
target. The optimal control law at the target is simply (u^ - 1, = 1) 

for the normalized system. However, this law is only valid for that 
single point in state space which defines the target. The state 
(x^ = 1 + e, x^ = 1 + <5) will have an entirely different law. 

In reality, we will consider the target to be some small region, 
not an infinitesimally small point. As far as the simulations are con- 
cerned, this small target region has already been determined by the 
quantization used in the Dynamic Programming. Whenever both states are 
within Ax (the quantization size) of the target, the simulator will 
begin to interpolate on the control law (u^ = 1, Ug - 1) for the target. 
This will cause the states to be slowly eased towards the target, which 
might be considered unacceptable, depending on how large a region of 
state space is involved. For this reason, these simulations determine 
approximate times at which the target point is reached, and exact times 
at which the target area is reached. 

Since all the Dynamic Programming solutions involve normalized 
values, both simulators must convert to the unnormalized system. All 
plotted figures in this chapter are normalized, for purposes of 
comparison. 

Although the specified goal is to talce a low- thrust starting state 

to the high-thrust target state in minimum time, the nature of an optimal 

* + 

control study (and in particular, the Dynamic Programming method) is 
that it is more concerned with state space and the time domain. As such, 


no plots of thrust are presented, although thrust plots would be very 
necessary it, frequency domain transfer function studies of linear sys- 
tems. It is assumed here that satisfactory output responses are 
obtained when the system constraints are not violated. 

The choice of the initial state for the following simulations is 
somewhat affected by the peculiarities of the DYNGEN (Model 1) computer 
program. Specification of the initial conditions is determined by an 
"off-design point" (see reference [10]), which is generally determined 
by specification of WFB and A8, Ideally, it would be desirable to 
choose a starting state as far from the target as possible, in order to 
demonstrate the usefulness of the global feedback control law. This 
often would involve specification of an extremely low WFB, and is easily 
accomplished. However, the control law at this point will be a much 
higher WFB , and results in convergence problems once the transient simu- 
lation has begun. The success or failure of the simulation to converge 
is strongly controlled by the TOLALL and DT variables (see reference 
[10] and the program inputs given in the Appendix) and amounts to much 
trial-and-error technique. Even these variables cannot totally control 
the convergence difficulties, and further changes in the DYNGEN program 
itself are sometimes required, as explained in [3]. 

Due to the above considerations, a somewhat high initial condition 
is used in the simulations, the feeling being that lower initial condi- 
tions would require an unacceptable amount of tampering .with the tapes 
on which DYNGEN is stored. For similar reasons, no control laws which 
were obtained from unconstrained models are tested on DYNGEN. 
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6.2 Model 2 Simulation Utilizing Model 2L2 (Unconstrained) Control Law 
The Model 2 controllers are implemented using the program shown in 

the Appendix. Figures 6.1 through 6.4 show the results of a simulation 
utilizing the Model 2L2 (unconstrained) control law, which was presented 
in Section 5.4. The starting state chosen corresponds to a thrust of 
approximately 80% of design thrust. The affects of interpolation on 
the control law are readily seen in Figure 6.1, showing the inputs ver- 
sus time. After 0.08 seconds, the controls are slowly eased towards 
their normalized design values (WFB = 1, A8 = 1). Figure 6.2 reveals 
an overshoot of approximately 100% for compressor speed (N ) and 30% for 

fan speed (N ) > certainly not a desirable response for a jet engine, 

£ 

Unfortunately, elimination of the constraints in the determination of 
the control law has resulted in these undesir-able consequences, as shown 
in Figure 6,3. Turbine inlet temperature (T^) has skyrocketed to 550% 
of its design value in only 0.02 seconds, and the surge margins have 
also reached intolerable levels. The state space trajectory, as shown 
in Figure 6,4, agrees remarkably well with the optimal trajectories 
which were analytically determined and presented in Figure 5.9 of the 
previous chapter. The time it takes to reach the target is also in 
agreement with the cost results as presented in Figure 5.10, approxi- 
mately 0.12 seconds. 

6.3 Model 2 Simulation , Utilizing Model 2 (Unconstrained) Control Law 
Figures 6.5 through 6.8 represent the results when utilizing the 

Model 2 (unconstrained) control law, as determined in Section 5.7. The 
same initial state is used here as for the linear controller of 
Section 6.2, and with remarkably close results. The control laws are 
slightly different, but the resulting state and constraint variable 
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FIGURE 6,5* Model 2 Simulation Utilizing Model 2 (Unconstrained) Control Law; 
br/. (Quantization) = .02 
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trajectories are quite similar to the trajectories produced with the 
Model 2L2 control law. This is certainly striking evidence that Model 2 
is a nearly linear system. 

6.4 Model 2 Simulation , Utilizing Model 2 (Constrained) Control Law 

The effects of constraints on the control' law can be quite clearly 
demonstrated. Figures 6.9 through 6.12 represent the results of the 
Model 2 simulation, utilizing the Model 2 (constrained) control law, 
as developed in Section 5.8. The starting state corresponds to a thrust 
of approximately 74% of the design value. The controls are considerably 
smoothed out, and Figure 6.10 shows that the state— time trajectories 
proceed to the target much slower and less abruptly, than was the case 

in the previous two sections. However, the constraints are now at 
acceptable levels, as evidenced by Figure 6.11. Furthermore, the control 
is riding both the turbine temperature and compressor surge margin con- 
straints from the time = 0.02 seconds to time » 0.16 seconds. This 
agrees with the constraint analysis as shown in Figure 5.21 of the pre- 
vious chapter. Unfortunately, the constraint limits, as given in 
equations (3.3-8), (3.3-9), and (3.3-10) are slightly exceeded, even in 
this simulation. This is not entirely unexpected, when considering the 
rather important fact that the effects of three state variables are not 
seen in these results. Recall that the optimal control law was derived 
(see Section 5.7), out of necessity, by employing linear approximations 
for states (3), (4), and (5) of Model 2 (see Table 2.5). The simulation 
of the controller as presented in this section uses no such approxima- 
tion, and hence, some variation is expected. Furthermore, the choice 
of initial conditions for states (3), (4), and (5) introduces yet 
another consideration, and in fact, the particular choices for this 
simulation were somewhat arbitrary. Regardless of these slight 
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deficiencies, the cost is in good agreement with the Dynamic Programming 
results, 

6.5 Model 1 Simulation Utilizing Model 2 (Constrained) Control Law 

In view of all the evidence accumulated in Chapters 2 and 5 which 
show the differences between Model 1 and Model 2, non-optimal results 
are expected when applying on Model 1 a control law which was derived 
from Model 2. Indeed, this is clearly the case, as demonstrated by 
Figures 6.13 through 6-16. The starting state corresponds to an off- 
design point of (WFB = 2.2, A8 - 2.95) on the DYNGEN simulator, and is 
the same starting point as was used in Section 6.4. While the constraint 
variables are within acceptable limits, the state-time trajectories 
resemble very slow ramp functions. It takes 0.34 seconds to reach the 
target area and, while Model 1 is known to react more slowly than 
Model 2, it is not expected that the cost be that high. Clearly the 
control law is not satisfactory . for use on Model 1. 

6.6 Model 1 Simulation Utilizing Model 1L2 (Constrained) Control Law 
Application of the Model 1L2 (constrained) control law produces the 

best results for Model 1. This is established by Figures 6.17 through 
6.20, using the same starting state as the previous two simulations 
(approximately 74% of design thrust). After 0.23 seconds, both rotor 

m 

speeds are within 1.0% of their respective design values, a significantly 
better performance than is provided by the Model 2 control law. It is 
somewhat slower, however, than the cost predicted by the Dynamic 
Programming results of Section 5.5 (.205 seconds). This is not dis- 
turbing, and perhaps quite satisfactory, when considering that Model 1 
is a 16th order nonlinear simulation. It must be expected that the use 
of a second order linear approximation in obtaining a control law cannot 
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possibly result in the exact prediction of the minimum time it takes the 
trajectory to reach the target. The control quite clearly rides the 
turbine inlet temperature and compressor surge margin constraints over 
most of the trajectory, in agreement with the constraint analysis of 
Figure 5.18. 


• CHAPTER VII 


SUMMARY 

The goal of this work was to obtain a global nonlinear optimal 
control for a two spool turbofan jet engine. Various models were 
developed, pursuant to this goal. Most important of these models was 
the nonlinear analytically-expressed Model 2, which correctly models 
most of the qualitative behavior of the jet engine, but which fails to 
achieve strong numerical agreement with the non-ana lytical Model 1 
simulator. The time optimal control program was then expressed in 
detail, and various constraints were added to the problem. Dynamic 
Programming theory and the Successive Approximations technique were 
explored, and applied to the problem of interest, while several improve- 
ments in the numerical programming were introduced. Analytical and 
numerical results were obtained for several models, both constrained 
and unconstrained. Finally, these results were tested on the two 
principal simulators. Model 1 and Model 2. 

Indeed, this study has successfully achieved time optimal feedback 
control laws for various models of the two-spool turbofan jet engine. 
Furthermore, valuable insight into the nature of the problem has been 
obtained, and much useful computer software has been developed. Unfor- 
tunately, all enthusiasm for the results achieved in this study must be 
tempered by the realisation that an optimal control law obtained from 
any model can only be as good as the model itself. For this reason, 
more work is needed to develop a better nonlinear analytical model, 
similar to Model 2 as presented in this study. 

As the accuracy of these models is further improved, more consid- 
eration should also be given to the details which so greatly influence 
the time optimal feedback control law: the determination of the 
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allowable controls, U; and the limits placed on the selected constraint 
variables. As the entire analytical problem formulation (model, con- 
straints, etc.) becomes closer and closer to the actual physical problem, 
more detailed solutions can then be obtained in the numerical analysis. 

In conclusion, this study should be viewed as one more step in the 
efforts to ac. :ve global optimal control laws for two-spool turbofan 
jet engines. It has accomplished much of its original goal, but leaves 
much more work remaining. 
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ANALYTICAL MODELS OF A JET ENGINE* , 

R. J. Leake and J. G. Comiskey 

Department of Electrical Engineering 
University of Notre Dame 
Notre Dame, Indiana 46556 

ABSTRACT 

This paper describes an algorithm for obtaining a nonlinear analytical modal of a 
jet engine from measurements of two equilibrium point values and the linearized A 
and B matrices at those points. The method is compared with more conventional proce- 
dures of interconnecting individual component approximations. 


"j 1. INTRODUCTION 

In this work we continue the study of nonlinear analytical models for a two spool 
^ 'turbo fan jet engine first reported in [1]. The model given in [1] is refined here 
| and compared with an entirely new model which is the main subject of this paper. In 
order to distinguish the various models, we make the following designations. 

™ Model 1. -A large, flexible generalized engine simulator called DYNGEN which has 
‘ ‘ been developed at NASA Lewis Research Center [2,3] and coded for a parti- 

cular hypothetical two spool turbofan engine. 

• 

” Model 2. An analytically expressed set of 5 nonlinear differential equations plus 

about 20 nonlinear static equations approximating the relationship between 
various engine variables. 

* 

• r Model 3. A relatively simple two-input, five-state model which can be generated 
automatically for any engine from data of two equilibrium points plus A 
• ■ • and B matrices by the algorithm to be presented below. 


;i 


2. DESCRIPTION OF MODEL 2 


A refined and updated version of the two-input, five-state, 
T analytical model presented in [1] is given in this section, 
-i tions are: 


I. 

I 


1 

I 


*= fuel flow (WFB) 

U£ * nozzle area (Ag) 

X^ *= compressor rotor speed (N^) 

Xj ** fan rotor speed (N^,) 

Xg - burner exit pressure (1^) 

X^ c after burner exit pressure (F^) 

Xg ** high inlet energy (U^) 

“ thrust (FG) 

Yg ** high turbine inlet temperature (T^) 
Ehe system is completely specified as follows. 


two-output nonlinear 
The variable designa- 


Constants 


J “ AJ ~ 778.26 
G " 32.174049 
R»M o .0252 
Y* » 1.4 
P 2 = 518.668 
I c “ PMIHP «* 3.8 
Ip «= PMILP = 4.5 


'COMB 


1.65 


V AFBN 49,77 
CVMNOZ = .9494 

CAPSF = 2116.217 

N_ DESIGN = XNHPDS 

N_ DESIGN = XNLPDS 

F 


10070 

9651 


H ** 20.71175 


C FC “ 

CM 

• 

C PF " 

.24 

C VB 

.20279 

C P1IT " 

.22589 

C PLT = 

.27938 


“**<*■%£ 


$ “ PCBLC = * 16 
a «= PCBLDU = .208 
8 - PCBLHP = .726 
Y « PCBLLP = .066 


Design Equilibrium Point 
WFB = 2.75 
A g » 2.9482558 


CNF » 1.02310 
T21 » 742.957 
CNC = .98730 
T 3 b 1467.47 
T 50 « 2103.47 
T 55 « 1789.15 
ZC » .81430 
SFC » .737071 


(Sea Level Static) 

N c - 11899.1 

N p *= 9873.95 

P, » 23.9299 
4 

U 4 « 586.467 
P ? *= 2.55142 

P 3 » 25.3522 
P 21 » 2.9960 
WFMAX = 203.123 
PCMAX = 10.270 
WG7 * 224.323 
WAC « 137.649 
ZF - .8333 
BYPASS = .609694 


FG = 13431.02 

T. «= 2892.04 
4 


WG50 = 134.364 
WG4 = 118.375 
WG55 - 135.818 
PFMAX - 3.3624 
WAF - 221.573 
WCMAX » 54.4151 
AWMAX = 1.5805 
WG24 « 88.5047 


State Equations 
30 7 

( 1 ) — - = (— ) 
v ; dt v ix 1 


( 2 ) 


f!i 

dt 


30 2 

<f> 




Ry* 

V COMB 

ry*t ? 

V 

AFBN 


[C PC HACf *'a-V + C ]'HT HC50<T 4- T 50>' 
Tfo tC pF WAF(T 2 -T n ) + C pLT W C 55(I 50 -T55'l 

IT.WA3 + WFB - WG4] 

(WG4 - WFB - WA3] 


(5) 


dU, 

<i 

dt 




COMU 4 


~~ [t 4 (wga - writ „ WA3 > + y*< t 3 wa3 " t 4 Wg4 * V 1+n)WFB}] 


Nonlinear Functions Required fm' Eouatlons 


(1) CNF - 


N. 


N„DHS T .GN 
F 


9651 


(2) T 21 a T 2 + 214.2732 CNF 2 - 4Q(A - 2.948255) 


(3) CNC = 




K„ 


n c design/i 21 7t 2 


10070/t 21 /51S.668 


(4) T 3 « T n + 743.2722 CNC 2 - 68 (A g -2. 948255) 

(5) 3^ ~ U ,?/ C VB 

(6) T 5q .= .727 T 4 

(7) P 3 • 1.05944 

(8) P 21 = -6.20568 4- .0129774 - .0185376 ? 3 

(9) . W RIAX 4- 261.01 CNF - 63.196 


(10) P yMAy = 3.516739 CNF - .23561 

r -2.31326SCP 
.502 1-e 


(11) WAF ^ 


FHAX“ P 21 


’] 


188.113 CRC' 


(12) W = 137.54 - 457.987 CNC + 564.325 CNC^ 

(13) AW Q = 6.492 - 4.9749 CNC 

(14) P ctJAS * 26.43184 - 89.0484 CNC + 109.7243 CNC 2 ~ 36.5756 CNC 3 


(15) VAC = 


'21 


-. 3662 (P cmax } 


/r 21 /51S."66S 


W + AW CMAX a ~ e 


(16) WA3 = Cl-*)VAC » .84 WAC 

(17) WG50 = 301.957 F4//r^ 

(18) WG4 = WGSO - 8$WAC - KG50 - .11616 WAC 

(19) WG55 = WG50 + y$WAC « WG50 + .01056 WAC 


(20) T k , - 106.002 - .86154 T En - .10458 CNC> ; T n ,f Ert 

2D t 21 50 

(21) T y *= .49661 T 35 + 205. 8S6 Py 
1121.786 PyAg 


21 


1 


(22) WG7 = 


s 


Nonlinear Equations for Outputs 


(1) FG - .02951 WG7/1934.415 T ? + 68553.365 + 2116.217 A g (. 53978 P ? -l) 


(2) 2C 


(3) 2F 


<y p 2i> - x 

p - i 
CMAX 

P 21 ~ 1 

P - 1 
r FMAX 


3. DIRECT METHOD 

We now consider a direct computer method for obtaining nonlinear models. Let 

X ■» f(x,u) - (1) 

with x an n vector and u an m vector denoting a dynamical system such as a jet 
engine, in which the state variables and parameters u remain positive throughout 
the system operation and there is a function g(u) such that for each equilibrium 
point 

f (x,u) = 0 •* + x » g(u) (2) 

The steady state system analysis involves the study of the function g(u). 

We propose to approximate the system (1) by 

£ A(x) [x - g (u) ] ' ^ (3) 

where A(x) is a square matrix which varies as a function of x. Notice that if x^ 
is an equilibrium point of (1), Xp = g(u Q ), then a linearization about this equili- 
brium point results in the linear system^ 

fix - Apdx + B D Su ^ ’ (4) 

and a linearization of the approximating system (2) at x^ - gCe^) results in 

dx *= AjCx^fix + [— AC^) (u D )36u * (5) 

Hence, the linearization of (2) will match the linearization of (1) if and only if 

a(x d )=a d> ‘ _a d 15 C V " b d ' . < 6 > 

Also, if A is invertible, as is often the case for jet engine models, equation (6) 
yields 

"Su B D ■ ^ 

The basic idea of the proposed direct method is to use the above developments toge- 
ther with function approximations for A(x) and g(u) to arrive at a nonlinear model. 
Since equilibrium points and linearized models at those points can be obtained by 
known algorithms, we shall use this fact. Our initial approach is to use just two 
equilibrium points,. say and The input information is thus 

W V V V .V V B w ^ 

r . 

where x~ and x‘ are design and off-design equilibrium points. A^ and A^ are the 
system A-matrices at these points and u~, u„, B_, B , are the associated parameter 
and Input matrices. Mathematically, w w 

tffCXjjiUjj) 3f ( x D» U D^ 


V W’ D* W* 


( 9 ) 


We shall employ a linear approximation for A(x), given by 


V x j 


*3 ~ X 


A(x) = A., diag - — + A^ diag - - (10) 

W \ V 11 "U. X 

J J J J 1 

in which diag (*) is a diagonal matrix which causes the jth column of A(x) to be 
interpolated linearly between the jth columns of A., and A^ with x. as the interpola- 
tion variable. • ** 

The parameter vector u is presumed to be made up of physical control variables, and 
parameters such as altitude and Mach number. The equilibrium function is to be 
approximated in a manner such that both the equilibrium values and the lineariza- 
tions of the approximating system (3) match those of system (1) at both x_ and Xy. 
This requires then that 

= *3) 8< V = *W (11) 

and also 


(U D ) = 


3u ' 


The method we propose here is to approximate each scalar component g.(u) of g(u) by 


a linear affine nower law form 


C nrl-1 C m+2 


£ l (n) = c l“l + - • - + c „"a + c 2mU“l “2 
for which the j tli partial derivative is 


nrhn , 

U m * C 2m+2 


3u^ C j + ^2mfl C nrf-j 


.. m-M rrfm 
u . . . u 
1 m 


■ 'Wow, if the variables are normalized and scaled such that 
• Ujj •= (1,1,... ,1) = JL Uy = (a,a, . . . ,a) = a_ 

then, the conditions of (11) and (12) can be put in the form 
. 3g i 

k « ~ (1) “ C, + c c 

j j 2nrfl nrfj ^ . 

3g. 2c. , " 1 

V *= - — ~ ( a ) = c + r c a 
srfj 3Uj j 2m*H m rfj 


S i a) Sc j + C 2nrH + C 2nrf2 


g i (a) aSc^ + c 2m+1 a 


and summing the first two of these over j yields 

Ik. c Ec. + c„ Ec , . 
i j 2nr)-l mf j 

2k.. *= 2c, + c„ . a J Ec ,, 
mfj j 2m+l m+j 


k 2m+l * IC j + C 2mfl + C 2nH-2 


/ 


I 

1 

J 


rc 


ale. + c_ ,_a 
j 2m+l 


2bH-2 

which is of the form 


m+j 


+ c 


2m+2 


*1 “ r l + r 3 r 2 


V 1 


5 2 " r l + V2 a 

5 3 " r l + r 3 + r 4 


( 18 ) 


e. « ar- + r_a + r, 

4 13 4 

which, incidentally, is the m=l condition also. This set of transcendental equa- 
tions is solved numerically for r jL’ r 2 ,r 3 ,r 4 anc * (16) Is then used to solve for 

each c j . In the event that (18) has no solution, a best fit is made on the second 
equation by varying r£ while the other conditions are satisfied exactly. 

4. ALGORITHM OF THE DIRECT METHOD 

In this section, we present an algorithm which serves to automate the process of 
finding a nonlinear model for a system 

x ** f (x»u) ( 1 ) 

to be approximated from by a normalized system. The algo- 

rithm will automatically perform the normalization and, hence, actually approximate 
the sjstem . . 

X *= f (x,u) » 


■where x, - 


i ~ X i^ X D * “ U j^ U D * approximating system is of the form 

i i 


' x ** A(x) [x-g(u)J 
where 


x„ - x. 
D i x 


- V 


4; . 

I" 

t; 

I 

*i r* 

4 ** ‘ 

I- 

| 
^ t 
< L*j * ; 


• A(x) = A^ diag ■ \ + A_ diag ¥ ^ 

i i i i 

and ■* 


L ** Ec*u* + c* ir «*** + c* 
x . j j 2mfl j J 2mF2 


fthere u* = a.u. + B-. 

j j j 3 


Algorithm I . 

i* In P“ C! yvVV n ' n . a . ‘'VvV\ 

2. Calculate: 


( 2 ) 

(3) 

(4) 

(5) 


Ajj » diagCl/x^) diag (x Q ) 

- diagCl/Xjj ) A^ diag(x D ) 
i 1 

Bp ** ‘JlagCl/Xp ) B d diagCu^ ) 

X i 


3. Calculate 


B « diagCl/Xj. ) B diag(u ) ' 

* •** * ^ORIGINAL PAGE id 

He ©E EOT QUALITY 

“3 " M Wi’ 

v^V ( VV 


4. Calculate: 

.i 


• k j 11 B D^lj k 2nrfl 1 ^ 

’ ^Vii fc L-2 " W " V 


± i 


5. Calculate 


j r i»- 


± *= l, . 


£- IkJ 

A j=i J 


i _ .i 
S 2 "" k 2nri*X 


« y fci, . 

2 

i ,i 

s 4 ° k 2nr« 


6. Go to Algorithm II. 

♦ 0 .. i i i i 

Ssndw s 2 7 S 3 s ^4 3 ^ 

„ i i i i 

Eecexve: r^, r 2 , r 3 , r^,Y 


7. Calculate: 
i 


C 2m4*l " r 3 


. i 1 

C 2nrf2 r 4 


*w - k i + T 


8. Output 


V***’ e 2mf2 

V 

A A 

V X 


V K 


1 . i i i 

c. = k. - r„ c , . 

2 J 3 mt: 


i *= 1, . . 


1 ** 1. . . 


.1 « !.■ 


. i - 1»< 


j 35 1 j a * • jin 


d. 38 X$ # • . >xi 


l- Algorithm II 

5^; X* input; ^X* ^2* ^3 * ^4 * ^ ^ 


2. Calculate: 


S 3" S 4 

*X “ S 1 " TT 

_ ®3“ s 4 

*2 " S 2 TT 




3. Minimize by line search: 


c 

r •'* 

[■■ ■ 


x-1 

X - 

a 

■ ? 

L 

r 

p 2 ' 

ax * 

’ p i 

l ~ a-l 
X - 


E 


x - 

a “1 
a-l 


for -10 x _< 10 , 
4. Calculate: 

r 2 3 


x ^ 0, 


r 2 

a -1 


x ^ 1 


Q37A£^ 


S A * “3^ a ~ a ^ 


- as 3 - r 3 (a -a) 

X 2 " " '" 1-a 

* 

5. Return to Algorithm 1.6 

' .5. NUMERICAL RESULTS 


* Ca l ~ S 2 + r 2 r 3 Ca ’ 2 - U) 


The algorithm of the previous section was applied to data obtained using DYNGEN with 
Xp ana u^ specified as in Section 2. An off— design point was obtained using u^ = 

(.72727, .72727), with the resulting normalized state iL = (.9000, .7897, .7381,. 
.9401, .9454). The normalized A and B matrices are 


.-3.8 -1.277 2.067, 

2.748 -5.39 1.585 

“ 377.9 49.51 -264.9 

31.26 139.39 -6.269 

v -176.5 23.91 -10.27 


-1.277 


A« 377.9 
V 31.26 


2.748 -5.39 

377.9 49.51 


-1.152 

-1.991 


1.448 

1.071 


86.807 78.91 

-88.69 27.83 


-37.4 


-246.7 


I.,V 


-4.744 -1.3888 3.2468 -1.4591 1.1969 

.82186 -26.726 2.5585 -1.8609 .45548 

475.73 137.55 -328.91 27.791 91.495 

-50.103 110.91 63.188 -116.69 8.2883 

-186.77 -67.682 -41.681 24.586 -243.23 


-.00259 .3553 

.2116 -.31618 

B - 12.54 -13.774 (1) 

W -.6201 -99.3 

157.78 6.84 

-.04546 .0013 " 

.0086 -.0121 
B = 2.434 -.613 (2) 

w .67865 -97.467 

203.44 .64755 


I _ Using the^parameter value a - .7, the c. coefficient which specify the equilibrium 
... function g(u) as in Section 4^ are given^by the matrix . .t 

T ? ~ p, 24267 -.00218 1.90082 8.09916 .02864 .73088 1 ” 

‘4** ' 1.01593 .85407 .89872 .66919 -.81879 -.05121 

-■* C *= .73445 ,10133 6.90586 3.09409 .011495 .15272 (3) 

.77234 -.35905 2.45867 2.87415 -07519S .66191 

.39503 -.27262 -3,44682 13.4468 .01838 .85921 

I- J 

’Tf, This matrix together with the values a =1.1 and B= 0.1 and the matrices tL and A^ 

Jm *: • * completely specify Model 3A. 

■ Another model which we will call Model 3B is easily obtained by using a linear 

Jjt affine approximation to g(u) such that gC^) = x^, g(u w ) = x. rf . Model 3B is specified 


-l 

a «* c , a 

2,31778, 8 

« -1.31778 and 

the co cffi 

.cient 

1 

".1553 

.0028 

1.0 

1.0 

0.* 

.8418 


.1619 

.1707 

1.0 

.1.0 

0. 

.6674 

C 113 

.5351 

-.1208 

1.0 

1.0 

0. 

.5857 


.5878 

-.49313 

1.0 

1.0 

0. 

.9053 


..2962 

-.2D99 

1.0 

1.0 

0. 

.9137 


. ( 4 ) 


I In o>der to compare the four models, a test point x^ far ^removed from x^ was chosen 
by Setting u„ - (.8,1), and calculating the equilibrium x . A step change to u = 
(1,1) then causes an acceleration transient back, to x^. The results of this com- 
T parlson for the rotor speeds are shown in Tabi§ 1® More detailed information is 
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f. A numerical algorithm for obtaining nonlinear analytical models for jet engines is 
.. presented,. The method is to separate th.e transient A(x) and equilibrium g(u) parts 
of thq system dynamics and approximate these using easily accessible data. The com- 
fll ponents of g(u) are approximated by a. linear affine power law form. The principal 
|. ' numerical difficulty is that all boundary conditions may be impossible to meet. The 
*•' algorithm then satisfies all but the second of the equations (18) in Section 3 and 
a ' fits the second as closely as possible. The variable Y is zero when an exact fit 
r occursj but, otherwise, it eauses a least squares fit on the derivative conditions at x^ 
X Of equations (16), Section 3, while matching the other conditions exactly. The ' 
jj free parameter in the linear affine approximation of Model 3B are handled similarly. 
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Abstract 


Recently, the CARDIAD (Complex Acceptability 
Region for DT Agonal Dominance) plot has been in- 
troduced and applied to the problem of designing 
dynamical prccompensat ion to achieve column dom- 
inance. This paper illustrates several basic 
features of the method while using it to design a 
single , low-order dynamical compensator which a- 
chicves dominance at five operating points of a 
realistic two-spool turbofan digital simulation. 

. 

1. Introduction 

The CARDIAD (Complex Acceptability Region for 
DIA gonal dominance) plot is a graphical technique 
for choosing dynamical compensators to achieve 
diagonal row or column dominance, as define. * by 
Rosenbrock (1]. Without essential loss of general- 
ity, the compensator is assumed to have its di- 
agonal elements equal to unity, and a typical 
CARDIAD plot describes the acceptable range of the 
reel and imaginary parts of the off-diagonal el- 
ements such that dominance is achieved. The basic 
graphical building block is the circle. Each 
circle represents the acceptable, range at a spec- 
ific frequency. Solid circles are drawn if ac- 
ceptable real and imaginary pairs correspond te 
points inside the circle, and dashed circles are. 
drawn if acceptable pairs correspond to points out- 
side the circle. Plotted as a function of fre- 
quency, these circles describe the acceptable range 
of the compensator element in question, considered 
as a complex function of frequency. 

Recently, CARDIAD plots have been shown to be 
an effective design tool in dynamical prccom- 
pensation of multivariable plants to achieve dom- 
inance [2, 3, 4, 5]. This paper focuses upon an- 
other aspect of the CARDIAD plot, namely its 
ability to assist with che classification of var- 
ious operating points of a nonlinear system with 
regard to their dominance possibilities and to help 
with the design of compensators which achieve dom- 
inance of nul tlple operating points. 

*Thts *^>rk was supported in part by the National 
Acronantjlcs and Space Administration under Grant 
KSG 304S and in part by thr National Scirnce 
Foundation under Grant ENG 75-22322. 
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, 2 . Specific Assumptions , 

Plant models used to construct the plors in 
the sequel have been generated from the general 
purpose digital jet engine simulator DYNGEN [6j 
under a load which provides behavior similar to 
that of the F-100 two-spool turbofan engine at sea 
level static conditions. The models have two in- 
puts, five states, and two outputs.. They are lin- 
earizations of DYNGEN obtained with the aid of the 
DYGAF.CD package [7] under development at NASA Lewis 
Research Center. Fhysical description of the 
states .can be found in the references [3]. The in- 
puts are fuel flow and exhaust area; the outputs 
are thrust and high turbine inlet temperature. 
Paraneteriza tion is accomplished through the nom- 
inal value of thr fuel flow WF3, which takes the 
five values 2.145, 2.31, 2.475, 2.64, and 2.75 
LBM/SEC, ranging from a low thrust condition to 
high thrust without augmentation. All the models 
have been normalized. 

Thus the plant transfer function matrix has 
two rows and two columns, and exhibits transfer 
functions of degree five in both numerator and 
denominator. Space limitations preclude their 
presence in this manuscript. 

Denote the plant by G(s). Then the issue is 
to select a precompensator K(s) in such a way that 
C(s) K(s) is column dominant [1]. In particular, 
it is desired to select £n£ K(s) so that column 
dominance is maintained over all five nominal fuel 
flow conditions. 


3. General CARDIAD Features 

If the origin of the CARDIAD plot for a given 
column is included by all solid circles and ex- 
cluded by all dashed circles, that column of the 
system is dominant without further compensation, 
in as much as the origin represents unity com- 
pensation. Thus, the eventual goal of compensation 
using the CARDIAD plot method is to arrive at a 
system where all the CARDIAD plots have this 
feature. If there exists a point on the real axis 
such that the point is included by all solid 
circles and excluded by all dashed circles in the 
CARDIAD plot for a given column, then the choice 
of the value of this point in the off-diagonal 



A* | *« VfT^AI ;•! r « •%<- t»" I ; f%£! : “Tf D 



I 



L 





ORIGIN A^- 

qE POOR- 


page is 
quality 


entry which the CARDIAD plot represents will make 
the column dominant at all frequencies. If there 
Is ;w> such point the CARDIAD plot describes the 
range of a frequency dependent off-diagonal entry 
which will nuke the column dominant. 


CARDIAD plots for two input, two output sys- 
tems have sore Interesting features. A circle at 
a specific frequency in the CARDIAD plot for one 
column will be solid if and only if the other 
column is dominant at that frequency. Thus, when 
a system is dominant at all frequencies, all the 
circles in the CARDIAD plots will be solid and all 
will contain the origin. Another interesting fea- 
ture is the effect of a column switch, that is, 
soil tiplicat ion by a matrix with the only non zero 
elements being ones on the off-diagonal. The 
effects of such a switching of the inputs are that 
all the solid circles become dashed, all the dashed 
circles become solid, and the shapes of the column 
one and two plots are switched. This fact will be 
used in the next section to achieve dominance in 
the various set point models. 


4. Design Example 


vut v v : r. 
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nominal value of the fuel flow increases. Since 
there is more margin for error in the lover nom- 
inal value of fuel flow models, a compensator 
which Is fit to a rough average of the five plots 
and which tends to be closer to the higher nominal 
value of fuel flow models, might achieve dominance 
in all five models. * 

I 

The average value of the center of the lowest 
frequency circle of the five plots is -9.81. This 
suggests that designing a compensator to fit the 
nominal fuel flow of 2.75 model which has as the 
center of the lowest frequency circle the value of 
--9.59 might achieve dominance in all of the mod- 
els. The second order funtion that was chosen is 

-.747s - 9 .59 
.014s* -.998s + 1. 

and the next compensator, K^fs), is 
1 --742s - 9.59 

K(s) - .014s 2 -.998s +1 

0 1 

4 

Thus, the overall compensation is K(s) given by 


The CARDIAD plots of the five uncompensated 
models are all very similar in shape. This great 
similarity suggests that one compensator might be 
found that will make all of the models column dom- 
inant. The uncompensated plots also show ,.iat a 
column switch would make the first column of each 
of the models dominant at all frequencies without 
further compensation. Thus, was chosen cO be 

0 1 

1 0 • 

** 

Fig ires 1-10 are the CARDIAD plots of G(s)K} for 
the five models. The repetition of the general 
shapes of the plots, which is unaffected by the 
column switch, is very apparent. The plots also 
show that the first column of each of the models 
is now dominant. This can be ascertained either 
by the fact that the origins of the column one 
plots are included by all solid circles and ex- 
cluded by all dashed circles or by noting that all 
of the circles in the column two plots are solid. 

To achieve dominance in the second columns of 
the models, it is clear that some sort of frequency 
depentent compensation will be necessary because 
there exist uo points on the real axes of the plots 
which lie inside all of the solid circles. A first 
choice of a function to fit the paths of the cir- 
cles could be a simple first order function which 
traces a semicircle through the complex plane as 
the frequency varies. However, it is desired that 
one such function be found that will work on all 
five of the models; so, a second order compensator 
will be used to fit better the shape of the circles 
at the hihger frequencies. Two things that should 
be noted about the shapes of the circles in the 
column two plots are that the circles tend to be 
larger for the lower values of fuel flow and that 
in general, the center of the lowest frequency 
(largest) circles moves toward the origin as the 


K(s) * 


-.742s - 9.59 
.oi4s 2 -.998s + 1 


Figures 11 - 20 are the CARDIAD plots of 
G(s)K(s) for the five models. It is clear that 
they are all dominant at all frequencies since 
all of the circles are solid and all include the 
origin. Thus, one compensator has been found 
which will make all five of the models considered 
in this paper dominant. 


5* Conclusions 


Through the use of CARDIAD plots, it has been 
possible to achieve dominance over a range of 
operating points of a jet engine simulation. The 
compensator given above also achieves dominance at 
all but a very narrow range of frequencies in the 
model of another operating point. The results 
suggest two things. First, using the CARDIAD 
plots as a guide, it could be possible to design 
a compensator which varies with the nominal value 
of the fuel flow and achieves global dominance 
over a wide range of operating points. This is 
currently being studied. Second, the repetitive 
shape of the CARDIAD plots over the range of 
operating points suggests that the CARDIAD plot 
might be a useful tool in the classification of 
operating points with regard to interaction. Such 
a feature could be quite helpful in analysis of 
which models to use over flight envelopes varying 
from sea level to high altitude and fron low 
through high thrust. 
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entry which the CARDXAD plot represents will ntake 
the column dominant nt all frequencies. If there 
is no such point the CARDIAD plot describes the 
range of a frequency dependent oft-diagonal entry 
vhlch will snake the column dominant. 

CARDXAD plots for two input, two output sys- 
tems have some interesting features# A circle at 
A specific frequency in the CARD I AD plot for one 
column will be solid if and only if the ether 
column is dominant at that frequency* Thus, vhen 
a system is dominant at all frequencies, all the 
circles in the CARDIAD plots vill be - solid and all 
will contain the origin- Another interesting fea- 
ture is the effect of a column switch, that is, 
Bultiplication by a matrix with the only non zero 
elements being ones on the off-diagonal* The 
effects of such a switching of the inputs are that 
all the solid circles become dashed, all the dashed 
circles become solid, and the shapes of the column 
one and two plots are switched. This fact vill be 
used in the next section to achieve dominance in 
the various set point models. 


• A* Design Example 

The CARDXAD plots of the five uncompensated 
models are all very 'similar in shape. This great 
similarity suggests that one compensator might be 
found that will make all of the models column dom- 
inant. The uncompensated plots also show Luat a 
column switch would make the first column of each 
of the models dominant at all frequencies without 


.further compensation. 


Thus, was chosen lo be 


0 

1 


1 

0 . 


Tigures 1—10 are the CARDXAD plots of G(s)Ki for 
the five models. The repetition of the general 
shapes of the plots, which is unaffected by the 
column switch, is very apparent- The plots also 
show that the first column of each of the models 
is now dominant* This can be ascertained cither 
by the fact that the origins of the column one 
*’ plots are included by all solid circles and ex- 
cluded by all dashed circles or by noting that all 
of the circles in. the column two plots are solid. 

To achieve dominance in the second columns of 
the models, it is clear that some .sort of frequency 
cfepentent compensation will be necessary because 
there exist no points on the real axes of the plots 
which lie inside all of the solid circles. A first 
choice of a function to fit the paths of the cir- 
cles could be a simple first order function which 
traces a semicircle through t?te complex plane as 
the frequency varies. However, it: is desired that 
one such function he found that will work on all 
five of the models; sq, a second order compensator 
vill be used to fit better the shape of the circles 
Jt the hihger frequencies. Two things that should 
be noted about the shapes of the circles in the 
column two plots are that the circles tend to be 
larger for the lower values of fuel flow and that 
in general, the center of the lowest frequency 
(largest) circles moves toward the origin as the 
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nominal value of the fuel flow increases. Since 
there is more margin for error in the lover nom- 
inal value of fuel flew models, a compensator 
which is fit to a rough average of the five plots 
and which tends to be closer to the higher nominal 
value of fuel flow models, might achieve dominance 
in all five models. 

1 

The average value of the center of the lowest 
frequency circle of the five plots is -9.61* Tills 
suggests that designing a compensator to fit the 
nominal fuel flow of 2.75 model which has as the 
center of the lowest frequency circle the value of 
— 9*59 might achieve dominance in all of the mod- 
els. The second order funtion that was chosen is 

-*742s - 9.59 


.0l4s^ -*99Ss + 1. 


and the next compensator, K^Cs), is 
l ~.742s - 9-59 

K(s) - ,OWs i 2 -,998s +1 

0 1 


Thus, the overall compensation is K(s) given by 

0 1 


K(s) 


-,742s - 9*59 
,014s 2 -,998s -r 1. 


figures 11 - 20 are the CARDXAD plots of 
G(s)K{s) for the five models# ft is clear that 
they are all dominant at all frequencies since 
*11 of the circles are solid and all include the 
origin. Thus, one compensator has been found 
which will make all five of the models considered 
in this paper dominant* 


5. Conclusions , 

Through the use of CARDTAD plots, it has been 
possible to achieve dominance over a range of 
operating points of a jet engine simulation. The 
compensator given above also achieves dominance at 
all but a very narrow range of frequencies in the 
model of another operating point* The results 
suggest two tilings, first, using the CARBIAD 
plots as a guide, it could be possible to design 
a compensator which Varies with the nominal value 
of tlie fuel flow and achieves global dominance 
over a wide range of operating points. This is 
currently being studied* Second, the repetitive 
shape of the CARDXAD plots over the range of 
operating points suggests that the CARDXAD plot 
might be a useful tool in the classification of 
operating points with regard to interaction. Such 
a feature could be quite helpful in analysis of 
which models to use over flight envelopes varying 
from sea level to high altitude and from low 
through high thrust. . a:- " 
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Fig. 3. WFB-2.31, G(s)K_ , Column 1. 
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Fig. 8. WFB=2.64, C(s)K., Column 2 
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Fig. 7. VFB-2.64, C(s)K, f Column 1 
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Fig. 15. WFB*=2.475, C(s)K{s), Column 1 


Fig. 12. WFB=2.1A5, G(s)K(s), Column 2 


Fig. 13. VFB-2.31, C(s)K(s), Colugm 1 
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Fig. 17. WFB=2.64, G(s)K(s), Column 1 


Fig. 19. WFB-2.75, G(s)K(s) Column 1 


Fig. 18. WFB-2.64, C(s)K(s), Column 2 


Fig. 20. WFB-2.75, G(s)K(s), Column 2 
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Foreword 


From the outset, the use of a Theme Problem has posed certain challenges. Authors 
from academic backgrounds tend to be in need of highly detailed information about 
plant and specifications, while workers in ixidustry and laboratories must often be 
satisfied with indirect information and sometimes with none at all. We have tried 
to arrange a reasonable compromise somewhere on middle ground. Our decision to 
select a problem related to a realistic modern turbofan engine had special ramifi- 
cations of its own, not the least of which was the fact that certain types of addi- 
tional data were precluded for proprietary or other reasons. We believed all along 
that the advantages of data realism outweigh the disadvantages of incomplete infor- 
mation. 

The chronology of the Theme Problem begins in late summer, 1976, during discussions 
with J. L. Melsa. Subsequent contacts with several potential Forum participants 
led to the drafting of a Tentative Theme Problem Description, which' was sent out to 
various >workers for critique in early 1977. When evaluations were in hand, a Theme 
Problem Description was prepared on March 1, 1977 and became the working document 
for authors preparing papers for the meeting. Communications with several addition- 
al researchers established the need for minor modifications and clarifications, 
which were decided at a committee meeting held during the Joint Automatic Control 
Conference at San Francisco in June, 1977. These decisions formed the basis for an 
addition Theme Problem Memorandum mailed to all participants on July 18, 1977. 

All these adjustments are included in the Final Theme Problem Description, which is 
included here. 

Any clarity which may be present in this final problem description is due in large 
part to the valued advice of many colleagues, among whom I must especially mention 
R. L. DeHoff, R. D. Hackney, B. Lehtinen, W. C. Merrill, J. L. Peczkowski, C. A. • 
Skira, and H. A. Spang, III. Credit for any and all obscurities must, of necessity, 
accrue to the author. 


M. K. Sain 

Notre Dame, Indiana 

September, 1977 
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FINAL THEME PROBLEM DESCRIPTION 


. ABSTRACT 

xhe rtic-me design problem should serve the Forum goals in at least three ways. First; 
it nhould help to unify the presentations and, thus, make them more useful for group 
study after publication. Second, it should help to make the Forum relevant to the 
present-day design world by focusing upon a real system of considerable current in- 
. .tercet. Third, it should help to delineate the state of computational readiness of 
the various design viewpoints, and so help to point out where additional numerical 
researches would be useful. 

Caveat: It is important to recognize the generally positive intent involved with 

the use of this problem. It is not intended that the theme problem usage degenerate 
4- 11 to a computational contest. 

1. INTRODUCTION 

*• • t 

A very important developing area for linear multivariable control has arisen because 
of recent increases in the complexity of aircraft turbine engines. Engines in use 
today have, essentially, the one - control variable of fuel flow, though some make use 
of a variable nozzle area which is not unlike the iris diaphragm that controls 
aperture settings in a camera. Engines in the not-so-distant future can be expected 
to permit control of vanes in the stator portions of the various compressor stages. 
Further down the development line are engines with enough variable geometry to re- 
ceive the informal designation of "rubber engines" by research engineers in the 
Indus try . 

It is widely accepted that the older, workhorse, hydromechanical control methods are 
not equal to these new tasks and that they will, therefore, give way to electronic 
digital control. The entrance of the digital computer opens up vast numbers of new 
design possibilities, which are now beginning to receive increased attention in the 
industry. The central role played by the aircraft turbine engines in civil and 
military aviation makes clear the economic import of these trends. It would be hard 
to select a more timely theme design example for comparison of linear control alter- 
natives than the jet engine. 

In the United States, a joint study is now underway on the Pratt & Whitney F100-PW- 
100 afterburning turbofan, a low-bypass-ratio, twin-spool, axial-flow engine. Spon- 
sored by the Air Force and by the National Aeronautics and Space Administration, 
this study focuses on the linear quadratic regulator theory, applied at multiple 
operating points in the control regime. 

One effect of the theme usage of such a plant in the NEC Forum should be a broaden- 
ing of the design discussion to include other design viewpoints as well. 

2. PLANT 


The numerical model of the jet engine is supplied in (A,B,C,D) form on Attachment 1. 
For the A and C matrices, note that columns 9-16 are listed below columns 1-8. This 
model is for zero altitude and for a power lever angle (PLA) of 83 degrees, which is 
near maximum non-afterburning power. The motivation for choosing this operating 
point comes from the fact that every engine has to pass through this condition, as, 
fo:: example, on takeoff. Also supplied is a list of the input, state, and output 
variables associated with this model. These two pages are taken from the report 
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R. J. Miller and R. D. Hackney, "F100 
Multivariable Control System Models/ 
Design Criteria," Pratt and Whitney 
Aircraft Group, United Technologies 
Corporation, West Palm Beach, Florida, 
November 1976. 


Because a number of the techniques which will be discussed at the Forum have graphi- 
cal aspects, it is planned to facilitate the inclusion of curves in the publication 
by limiting the plant to three control inputs. In consultation with numbers of our 
Theme Problem Advisory Committee, we have selected U^, and as these inputs. 

Workers who feel an absolute necessity to use all five inputs are welcome to do so; 
however, we would ask that in such a case they provide a comparison of the effect 
of using five inputs over and above that of using only three. This request is 
designed to increase the comparability of the various design results. 


I 


* 


Actuator information for the three control inputs is given in Attachment 2. Also 
provided is information associated with the actuation of U,, if that input is used 

J 3 J • AJ J TT TT J IT T? • 1 1 1 J TT 1 3 • J J ' a VT TT 


in addition to U-, U^, and U 
.and a servo time constant of 


Finally, should U_ be used in addition to U- , « 9 , 

A A A 1 . 3 ^ J 1. • 


rate 


m , Cl 9C1 VU LXUiC CUUO CQU L UJL 0.02 sec. can be assumed for actuation, 
limits on the actuators can be noted, as in Table A. . 


* U 2 
Various 


Table A 

Actuator Rate Limits 


U 1 

15,800 

(lb/hr) /sec 

U 2 

3.6 

2 

Ft /sec. 

°3 

48 

Deg/sec. 

U 4 

40 

Deg/sec. 

The actuators have some limits 

, also. 

which will be 


be assumed that the limit is + 6°. On U ? , a limit arises because the nozzle area 


is pretty well down to its minimum at this operating point; the limit is assumed to 
be about 1 square feet in that direction. 


The Theme Problem models are in absolute, unnormalized form, without any mention of 
the set point values. This makes it difficult to size inputs. The committee worked 
out a proposal to supply "ballpark" set point values so that the model could be 
normalized. Unfortunately, it was not possible to obtain even such approximate in- 
formation. 


A consequence of this fact is that the absolute rate limits of Table A have meaning 
only in relationship to the size of reference commands assumed. Because we are un- 
able to supply the suggested reference command, the effect of actuator rate limits 
can be treated only hypothetically; and we have to leave the issue of whether to do 
this, and how to do this, in the hands of the authors. 


Turning now to the sensed variables, we have available X 


and (X 10 + 


»vw uu cut V ui .lciuo.u. 1 J wt nave a vaxxauxv '*1* ) <** ^ y auu 

X- ,) , the last of which is denoted FTIT for "fan turbine*" inlet temperature." 

O ~ i. J j • J ^ 1 J _ j 3 J _ rp_Vl - Tj 
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Sensor time constants in seconds are listed in Table 3 


1 




Table B 

Sensor Time Constants 


original page is 
OB EO OR QUALn* 

X 1 

X 2 

0.03 

0.05 


X 3 

0.05 


X 5 

0.05 


Sensing of FTIT is a bit more elaborate and is indicated on 


Attachment 


3. 


3. ENVIRONMENT 


Measurement noise is on the order of 1%; and state noise is negligible. Therefore 
it is not planned to supply any noise data. Authors wishing to make noise studies 
must make their own assumptions. This is not unrealistic for the present stage of 
discussion. Though some techniques may well make use of observers or dynamical out- 
put feedback* no formal stress on filters is anticipated. The Forum, then, is 
visualized primarily as a control meeting, although contributed papers in the stoc- 
hastic area will be accepted if they contribute to the Forum theme. 


Practice in the industry involves the use of multiple linear models at various oper- 
ating points from sea level to high altitude and from low to high thrust. As oper- 
ation transitions from the neighborhood of one operating point to the neighborhood 
of another, these models change in consonance with some physical variable. Para- 
meter Vdiiation is, therefore, an aspect of design. 

But publicly available neighboring linear models are not near enough to the Theme 
Problem model to provide meaningful data on parametric variation. This fact, com- 
bined with lack of set point information, led the committee to suggest a 5% change 
in eigenvalues as one, hopefully useful, measure of such variation. Because normal- 
ization of the model is a similarity transformation, this characterization is in- 
dependent of set point. 


4. REDUCED ORDER MODELS 

Approximate eigenvalues of the Theme Problem plant are -577, -176, -59.2, -50.7, 
-47.1, -33.7, -21.3 + i.822, -17.3 + i4.78, -19.0, -b.71 + 11.31, -2.62, -1.91, 
-.648. It is the nature of the jet engine control problem that these can usually 
be well identified with physical variables. For example, -.648 associates with X^g, 
-1.91 associates with , -2.62 associates with X 0 , and so forth. X^ is related 
to the eigenvalue pair -6.71 + il.31. This type of information can be deduced 
from a study of the eigenvectors corresponding to a particular eigenvalue. It can 
be expected that actuator modes, such as that involved with fuel flow, will enter 
into this list. Some discussion on his point can be found in R. L. Delloff and 
W. E. Hall, Jr., "Design of a Multivariable Controller for an Advanced Turbofan 
Engine, 11 Proceedings 1976 IEEE Conference on Decision and Control, page 1002. 

In the Interest of offering some assistance to authors who might be having computa- 
tional difficulty with the full size problem, the following reduced model has been 
made available by Dr. DeHoff of Systems Control, Inc. (Vt.). It is a model which 
neglects sensor dynamics, augments the plant by the dominant actuator dynamics, and 
then reduces to fifth order. The resulting five states are 

A 

= Fan Speed (rpm) 


5 


X£ " Compressor Speed (rpm) 

A 

- "Augmentor" Pressure (psia) 

X. ■= Fuel Flow (lb. /hr.) 

4 

A 

X^ = Burner Pressure (psia) 

A 

Note that the "Augmentor" Pressure is not to be identified with X^; the quanti- 
ties are not defined at the same physical location. Note also that X,. was not 
one of the original states. 

Remark : The U„ Actuator diagram shows a Servo System gain of 2.4. It has come to 

our attention chat a more realistic number for this gain would be about 12.0. The 
effect of this gain change is to take the dominant CIW position actuator eignevalue 
from a location of high dominance in the overall plant -actuator system to a location 
of considerably less dominance. It is not necessary for authors to make this change 
if they ( have already completed their calculations, inasmuch as the 2.4 gain apparent- 
ly is one of those "glitches" which crept in an uninvited manner. Some authors may 
choose to compare the effect of the gain 12.0 with the gain 2.4, if time and space 
permits. We have included this remark here so that the reduced order model, which 
has the same controls and outputs as the full size system, may be more understand- 
able. 


A (5 x 5) 


-.3245E+01 

- . 2158E+01 

-. 9155E+03 

.5731E+00 

.1342E+03 

.1 '422+01 

-. 5941E+01 

-. 2816E+03 

. 1897E+00 

•5705E+02 

.1685E-01 

-.2554E-01 

-.1003E+02 

.7994E-02 

•5807E+00 

.0000 

.0000 

.0000 

-.1000E+02 

.0000 

-.2163E+01 

.6862E+01 

. 7405E+03 

.1195E+01 

-.1715E+03 

B (5 x 5) 
.1432E-01 

-.'3553E+03 

-.9906E+02 

-.1549E+02 

. 2220E+05 

. 2871E+00 

. 7286E+03 

.2514E+02 

-. 6487E+02 

. S122E+04 

-.2469E-02 

-.1030E+03 

. 6333E+00 

-.3213E+00 

-.7418E+02 

. 1000E+02 

.0000 

.0000 

.0000 

.0000 

-.1311E+00 

.3295E+03 

-. 250QE+02 

•6257E+02 

-.6445E+05 

C (5 x 5) 

•1662E+01 

-.1768E+01 

. 7999E+02 

-.1890E+00 

.3771E+02 

. 1383E-01 

.3142E-05 

-. 1060E-01 

•1289E-03 

-.1839E-06 

.1694E+00 

-.1129E+00 

-.4959E+01 

.7386E-01 

-.1835E+00 

•7590E-04 

.3269E-05 

-.1477E-01 

.2284E-05 

.4315E-04 

-.4859E-04 

.1381E-03 

.1140E-01 

.1951E-04 

-.2688E-02 

D (5 x 5) 

•1302E+00 

.1992E+03 

.4802E+02 

-.1503E+02 

.1083E+05 

.1449E-06 

.3395E+00 

. 6806E+00 

.2812E-03 

.3204E-03 

.2967E-01 

.7927E+02 

.2567E+01 

-.7631E+00 

. 2066E+04 

.1046E-05 

-.7720E-02 

-.5814E-02 

.1157E-03 

.6605E-01 

-.8395E-05 

-.7897E-02 

-.6841E-03 

-.9643E-03 

-.2815E+00 
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5. SPECIFICATIONS 


The overall viewpoint of the controller is quite simple. The pilot has one lever, 
which we might intuitively call the throttle and which sets what is called in the 
industry the "power lever angle." Basically, the pilot increases the lever angle 
to obtain more thrust. All the ether variables must be controlled so as to achieve 
the new thrust quickly, but without overshoot and without violating some important 
physical considerations. An example of one of these is the temperature at the inlet 
to the "high" turbine just aft of the burner. This temperature is ordinarily 
scheduled very near its maximum safe value, and temperature increases are not wel- 
come because the turbine elements are thin, respond very fast, and can be permanent- 
ly damaged or create a need for more frequent engine overhauls. Another example of 
a constraint is the various undesirable stall conditions in the compressor. 

This problem comes down to us in the following form. Assuming a step change in 
power lever angle, we want to move the engine to a slightly different operating 
point in the above described acceptable dynamic fashion. The power lever angle 
change is converted by a master engine scheduler into a reference input for our 
linearized feedback model. The nature of this reference input is not highly specif- 
ic. Step inputs are commonly studied. It is not likely that highly detailed infor- 
mation about these references will be available, but we can try to firm up any par- 
ticular issues which may be crucial to one paper or another. The exact nature of 
these references gets one into the exact nature of the schedulers. It does not seem 
too productive in a linear meeting to go very far into such "global" issues. If 
greater reference variety is needed, it can probably be safely assumed. It would 
be good, however, if each paper tried to discuss at least the reference step. 

For purposes of design, we can group the variables into two families. Y^, Y^, , 

and X 2 are desired to respond fast without overshoot. Y^ should not decrease more 
than .05; Y^ should not decrease more than .15. 

Remark: The decrease limits on Y^ and Y are to be regarded in the same spirit as 

the actuator gain change in the preceding section. If calculations are complete, 
there is no requirement to incorporate it. Some authors may wish to study its 
effect, however. 

6. VIEWPOINT 

We believe that the theme problem should appear in each presentation as the major, 
and probably the only, illustration of the particular design methodology being 
described. We visualize each paper as an exposition of design viewpoint, with jet 
engine illustration. We do not visualize the paper as an exposition of jet engine 
design. In other words, the theme problem will be an apparent thread through the 
fabric of the Forum, but the pattern of the fabric will be set by the various linear 
control alternatives as entities in themselves. Put in yet another way, the Forum 
is on "Alternatives for Linear Multivariable Control" and is not upon "Various 
Approaches to Jet Engine Control." 
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F100 MODEL 

ALT-0.0 PLA=83 




THE A MATRIX 






-4.328 

.1714 

5.376 

401.6 

-724.6 

-1.933 1.020 

,-.9820 

-.4402 

-5.643 

127.5 

-233.5 

-434.3 

2o.59 2.040 

-2.592 

1.038 

6.073 

-165.0 

-4.483 

1049. 

-82.45 -5.314 

5.097 

.5304 

-.1086 

131.3 

-578.3 

102.0 

-9.240 -1.146 

-2.408 

.8476E-02 

-.1563E-01 

.5602E-01 

1.573 

-10.05 

.1952 -.8804E-02 

-.2110E-01 

.8350 

-. 1249E-01 

-.3567E-01 

-.6074 

3.7.65 

-19.79 -.1813 

-.2S62E-01 

.6768 

-.1264E-01 

-.9683E-01 

-.3567 

80.24 

-.8239E-01 -20.47 

-.3928E-01 

-.9696E-01 

.8666 

16.87 

1.051 

-102.3 

29.66 .5943 

-19.97 

-.8785E-02 

-. 1636E-01 

.1847 

.2169 

-8.420 

.7003 . 5666E-01 

6.623 

-.1298E-03 

-. 2430E-03 

. 2718E-02 

. 3214E-02 

-.1246 

. 1039E-01 . 8395E-03 

. 9812E-01 

-1.207 

-6.717 

26.26 

12.49 

-1269. 

103.0 7.480 

36.84 

-. 2730E-01 

-.4539 

-52.72 

198.8 

-28.09 

2.243 .1794 

9.750 

-.1206E-02 

-. 2017E-01 

-2.343 

8.835 

-1.248 

. 9975E-01 . 8059E-02 

.4333 

-.1613 

-.2469 

-24.05 

23.38 

146.3 

1.638 .1385 

4.486 

-.1244E-01 

.3020E-01 

-.1198 

-.4821E-01 

5,675 

-.4525 19.81 

.1249 

-1.653 

1.831 

-3.822 

113.4 

341.4 

-27.34 -2.040 

-.6166 

.9990 

1.521 

-4.062 

9.567 

10.08 

-.6017 -.1312 

•9602E-01 

11.32 

10.90 

-4.071 

-.5739E-01 

-.6063 

-.7488E-01 -.5936 

-. 9602E-01 

-.9389E-02 

.1352 

5.638 

.2246E-01 

.1797 

. 2407E-01 1.100 

. 2743E-01 

-3.081 

-4.529 

5.707 

-.2346 

-2.111 

-.2460 -.4686 

-.3223 

. 2090E-02 

-.5256E-01 

-.4077E-01 

-.9182E-02 

-. 8178E-01 

. 3428E-01 . 4995E-02 

-.1256E-01 

-.1953E-01 

-.1622 ‘ 

-.6439E-02 

-. 2346E-01 

-.2201 

-. 2514E-01 -.3749E-02 

-.3361E-01 

.1878E-01 

-.2129 

-.9337E-02 

-.3144E-01 

-.2919 

-.3370E-01 . 8873E-01 

-.4458E-01 

. 2253E-01 

.1791 

. 8371E-02 

. 2645E-01 

.2560 

. 2835E-01 -.3749E-01 

.3635E-01 

-49.99 

. 6760E-01 

39.46 

.4991E-02 

.8983E-01 

. 5349E-02 .0 

.1372E-01 

-.6666 

-.6657 

.5847 

. 6654E-04 

.1347E-02 

. 7131E-04 .0 

.2057E-03 

.2854 

2.332 

-47.65 

.3406 

3.065 

.3624 -.4343 

.4681 

-9.627 

-9.557 

38.48 

-50.01 

.1011 

. 1203E-01 -.4686E-01 

. 1715E-01 

-.4278 

-.4245 

1.710 

-2.000 

-1.996 

.5J49E-03 -.1999E-02 

. 7544E-03 

-4.414 

-4.354 

17.66 

-3; 113 

-3.018 

-19.77 -.4999E-01 

.1509E-01 

-.1127E-02 

-. 6760E-02 

. 1835E-01 

-.9981E-03 

-.1347E-01 

-.1070E-02 -20.00 

-.2057E-02 

.5004 

-.1437 

-2.416 

-.1073 

-1.078 

30.63 19.89 

-50.16 



ATTACHMENT 1 


THE B MATRIX 


.4570E-01 -451.6 

-105.8 

-1.506 

851.5 

.1114 -546.1 

-6.575 

-107.8 

3526. 

.2153 1362. 

13.46 

20.14 

-. 6777E+05 

.3262 208.0 

-2.888 

-1.653 

-269.1 

. 9948E-02 -98.39 

.5069 

-.1940 

-94.70 

. 2728E-01 71.62 

9.608 

-.3160 

-184.1 

.1716E-01 71.71 

8.571 

.7989 

-515.2 

. 7741E-01 -141.2 

-.8215 

39.74 

1376-. 

. 3855E-01 -7.710 

-.4371E-01 

-.1024 

-6684. 

. 5707E-03 -.1144 

-. 6359E-03 

-.1432E-02 

-99.02 

5.727 -1745. 

-8.940 

-17.96 

•8898E+05 

.1392 -24.30 

-.2736 

-.3403 

-6931. 

•6172E-02 -1.082 

-.1183E-01 

-.1452E-01 

-307.7 

.6777E-01 16.60 

.3980 

. 2311E-01 

-2588. 

. 1880E-02 9.147 

-.8241 

.8984E-01 

-32.31 

.1677 435.8 

-89.94 

4.900 

-295.5 


THE C MATRIX 


.4866 

-.6741 

5.392 

95.42 

•1383E-01 

. 2789E-05 

.0 

.0 

.0 

.0 

.0 

.0 

. 7418E-04 

•5496E-05 

.4790E-05 

.1478E-03 

. 1538E-04 

.1201E-03 

-. 2579E-02 

-. 1609E-03 

.5195 

.8437 

-1.863 

•5709E-01 

.0 

.0 

.0 

.0 

.0 

.0 

1.000 

.0 

.3434E-05 

. 2727E-04 

.1128E-05 

.4002E-05 

-.3732E-05 

-. 2996E-04 

-.1234E-05 

-.4380E-05 

THE D MATRIX 



-. 6777E-01 

-420.5 

32.97 

-1.824 

.1282E-03 

.3353 

. 6804 

-.5605E-04 

.0 

.0 

.0 

.0 

.1030E-05 

-.1193E-01 

-.5806E-02 

•6015E-04 

.8109E-05 

. 2328E-01 

.1178E-03 

-. 5538E-02 


24.03 10.52 .8190 -.4492 

-.1081E-01 -.5545E-04 .4722E-04 .0 

.0 .0 .0 .0 

-.1504E-01 -. 6503E-04 .8820E-04 .4999E-05 

. 1618E-01 -.1071E-02 -.9561E-04 -.5503E-05 

.4815 3.428 2.161 .7681E-01 

.0 .0 .0 .0 

.0 .0 .0 .0 

. 3673E-04 . 4290E-05 -.4958E-05 .5609E-05 

-.4024E-04 -.4721E-05 .5324E-05 -.6103E-05 


1245. 

-.1199E-01 

.0 

.4463E-01 

-.1039 


Engine State Variables 

X x «= Fan Speed, SNFAN (N^ - rpm 

X2 = Compressor Speed, SNCOM C^) - rpm 

X^ = Compressor Discharge Pressure, P - psia 

X. = Intnrturbine Volume Pressure, P . c - psia 

X c = Augmentor Pressure, P , - psia 

j \ /m 

Xg = Fan Inside Diameter Discharge Temperature, T^ - °R 

Xj = Duct Temperature, T^ “ ° R 

Xg = Compressor Discharge Temperature, T fc g - °R 

Xg = Burner Exit Fast Response Temperature, T^^^ - °R 

X^Q = Burner Exit Slow Response Temperature, T , “ ° R 

X.. = Burner Exit Total Temperature, T fc4 - °R 

X^2,*“ Ran Turbir a Inlet Fast Response Temperature, T fc4 - 

X^g = Fan Turbine- Inlet Slow Response Temperature, T^ c ^ o - 

X^ = Fan Turbine Exit Temperature, T^ “ ° R 

X^ = Duct Exit Temperature, T t g c - °R 

X^g = Duct Exit Temperature, T t 2 m ~ ° R 


°R 

°R 


Engine Inputs 

U 1 = Main Burner Fuel Flow, WFMB - lb/hr 

1 2 
U2 ® Nozzle Jet Area, A 4 - ft 

Ug = Inlet Guide Vane Position, CIW - deg 

= High Variable Stator Position, RCW - deg 

Ug = Customer Compressor Bleed Flow, BLC - % 


Engine Outputs 

Y^ » Engine Net Thrust Level, FN - lb 
Y2 = Total Engine Airflow, WFAN - lb/sec 
Yg = Turbine Inlet Temperature, T^ - °R 
Y 4 - Fan Stall Margin, SMAF 
Yg = Compressor Stall Margin, SMHC 
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ATTACHMENT 2 


Metering Valve 


Pump Controller 


Fuel Flow 
S* 

Request 


1 


1 

0.02S + 1 


0.1S + 1 


Fuel Flow 

-> 

To Engine 


Actuator 


Servo System 


0.01S + 1 


Air Motor Dynamics 

I 1 


JLs* 

25 

4_ e 

a> 2 
n 

i o 

0) 

n 

0) = 

6Hz 

n 


C = 

0.56 


U 2 Actuator 


Stepper Motor System 


0.02S + 1 


Servo 

System 

2.4 


0.01S + 1 


Power 

Cylinder 


40 


0.01S + 1 




nder 

RCW Position 
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ABSTRACT 

The determinant of return difference establishes a crucial link between open and 
closed loop characteristic polynomials in multivariable feedback control systems. 

As a result, Nyquist constructions on this determinant carry important design infor- 
mation. One way to extract this information is by achieving diagonal dominance. 

This paper presents a method which uses dynamical input compensation to achieve 
column dominance. Application to the Theme Problem is included. 
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INPUT COMPENSATION FOR DO* '...<CE OF TU1U10FAN MODELS 

R. M. Schafer and M. K. Sain 

Department of Electrical Engineering 
University of Notre Dame 
Notre Dame, Indiana 46556 
U.S.A. 

ABSTRACT 

The determinant of return difference establishes a crucial link’ between open and 
closed loop characteristic polynomials in multivariable feedback control systems. 

As a result, Nyquist constructions on this determinant carry important design infor- 
mation. One way to extract this information is by achieving diagonal dominance. 

This paper presents a method which uses dynamical input compensation to achieve 
column dominance. Application to the Theme Problem is included. 

1. INTRODUCTION 

Recent advances in the generalized Nyquist theory for linear multivariable feedback 
control systems have brought about very substantial new opportunities for research 
in the area of frequency domain control design. Most of these advances are predi- 
cated upon the relationship between closed loop and open loop characteristic poly- 
nomials — as embodied in the determinant of return difference. Features of the Nyquist 
diagram of this determinant are important aids to control system design. 

It is apparent that a diagonal return difference will decompose <-he return differ- 
ence determinant into a product of its diagonal elements, thus reducing a multi- 
variable problem to classical single- input, single-output form. Less apparent, but 
of much greater practical significance is the fact that an approximately diagonal 
return difference can have essentially the same reducing effect on a multivariable 
problem, when regarded from a generalized Nyquist viewpoint. The best known of 
these approximately diagonal conditions has come to be described as diagonal domi- 
nance. A productive design strategy can be mounted, therefore, in two steps. First, 
achieve diagonal dominance; second, apply classical single-input, single-output 
techniques [1]. 

Unfortunately , methods to attain diagonal dominance have been rather slow to advance. 
For the most part, they have been restricted to the selection of constant real com- 
pensators, the entries of which are typically obtained by procedures of optimization 
that do little to preserve some of the classical advantages, such as insight, affor- 
ded by the frequency domain approach. Much work needs yet to be done on the theory 
of attaining diagonal dominance by use of frequency dependent, dynamical compensation 

This paper considers the application to the Theme Problem of a useful new design aid 
called the CARDIAD Plot. In its present form, this method deals with the design of 
a dynamic precompensator for the plant, in such a way that column dominance is 
achieved. An important feature of the approach is the enhancement of designer in- 
* sight toward the coupling present in a plant. 

Section 2 introduces the CARDIAD method for two-input, two-output plants, and Sec- 
tion 3 provides an illustration of certain basic features of the method, in the con- 
text of a jet engine plant related to the Theme Problem. Section 4 gives a genera- 
lization of the idea to three inputs and three outputs, and Section 5 applies these 
*• results to the Theme Problem. Conclusions appear in Section 6. 
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2. GRAPHICAL APPROACH 


I 

I 
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The i tb column of a matrix Z(s) is said to be dominant if 

- j, bi (s) | * 0 

j^i 

for all s on a Nyquist contour D. A similar definition can be made for row domi- 
nance. 


( 1 ) 


I For a two- input, two-output system, Eq. (1) can be equivalently written 

1 . . |2 I ..I 2 A a f m 

if j 


i 

I 


z ±1 (s) 2 - | Zji (s)| 2 > 0 


(2) 


for all s on D. - 


Consider a two-input, two-output system having only precompensation. The open loop 
transfer function of the system is 


Q(s) = G(s)K(s) . 
<(s) 

K(s) 


_ Let K(s) be restricted to the form 

I r 1 


s 


o 2 (s) 


a^(s) 


(3) 

(A) 


' Since any matrix having nonzero entries on its main diagonal may be put into this 
form by multiplication with a diagonal matrix, and since, multiplication by a diago- 
nal matrix does not affect dominance , this can be done ‘without essential loss of gen- 
erality. 


(5) 


evaluated at 

a specific 

frequency m. Then 

v u + y 

r 12 + 1 12 3 ’ 


1 x 2 + y 2 3 " 

_ r 21 + ± 21 i 

r 22 + 1 22 3 _ 


X 1 + y l 3 1 


Q(jw) = 

Performing the indicated matrix multiplication, the four entries in the matrix 

Q(s)| are 
|s=> 

’ll “ r ll + £ 11 3 + (r 12 + l 12 3 ’ tx l + y l 3 > • . 

"12 " r 12 + *12 3 + <r U + hi 33 '’ 1 ! + y 2 3) ’ 

’21 “ r 21 + *21 3 + (r 22 + i 22 J>(x l + y l 3) > ■ 

• ’22 ’ r 22 + *22 3 + (r 21 + ^l 3 3 (x 2 + y 2 3, ‘ 


From Eq. (2), the first column of Q(s) will be dominant if 

J ,2 , ,2 


I 


*11 


- q 


21 


> 0. 


( 6 ) 

(7) 

( 8 ) 

(9) 

( 10 ) 


.Performing the indicated subtraction results in what will be referred to as the domi 
nance inequality for column 1. The form of this inequality is 

*l^ X l’ y i^ “ 3X 1 + ay l + 2bx l + 2cy l + d °» 

* where the constants arc defined as 
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2 A ,2 2 ,2 

a " r i2 + i 12 ~ r 22 “ *22 


( 11 ) 
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i ll “ r 21 


12 11 21 22 11 12 22 21’ 7"*''WAC 

2 2 2 2 0S fOOR OH. * « 

d « . .*2 2 2 ' ^^ALlTy (14) 

d r ii + hi r 2i hr • 1 r k ; 

Note that each constant is composed of complex field elements which come from evalu- 
ation of G(s) at a specific frequency w. 

The function f^x^y^ is a paraboloid in three-space and is normal to the 

plane. If this paraboloid intersects the x^-y^ plane, the intersection will be a 
circle. Standard maximum-minimum analysis gives that the maximum or minimum of the 
dominance function occurs at 


-b/a 


-c/a 


To determine if the point that was found is a minimum or a maximum, the hessian is 
formed. If the hessian is negative definite, the point found is a maximum. If the 
hessian is positive definite, the point found is a minimum. The hessian of the domi- 
nance equation for column one is 


2a 0* 

0 2a 


so that the second derivative test reduces to a test on the sign of a. 

Proceeding from this analysis, there are four possible cases. The point that was 
found vas a positive maximum, positive minimum, negative maximum, or negative mini- 
mum. The two cases that are of interest are the positive maximum and tne negative 
minimum since it has been shown [2] that the other two cases cannot occur. In each 
of the cases of interest, the positive maximum and the negative minimum, there is 
an intersection of the x -y^ plane. Recalling that the column will be dominant if 
^1^ X 1*^1^ "*' S P os * t ^- ve » t “ e analysis of the two cases is as follows. In the positive 
maximum case, the values of x^ and yj_ which will result in solution of the dominance 
inequality are those points which lie inside the intersection of f^Cx-py^) and the 
x l**^l P^ ane » that is the circle which is the solution of f^(x^,y^)=0. In the nega- 
tive minimum case, the choices of x-^ and y^ which result in solution of the domi- 
nance inequality are those points which lie outside the circle of intersection. 

Thus, the intersection of the dominance function f^(x^,y^) for column one and the 
Xj-y^ plane defines the acceptable range of x^ and y^ such that the system will be 

dominant in the first column at the specific frequency at which the analysis was per- 
formed. In like fashion, the second column of the system may be analyzed, and the 
acceptable choices of and y^ may be determined. 

If this dominance analysis is repeated over a range of frequencies, and the result- 
ing circles of intersection plotted, a CARDIAD (Complex Acceptability Region for 
PI A gonal Dominance) Plot is produced. A solid circle is drawn if the acceptable 
choice of x and y lie inside the circle, and a dashed circle is drawn if the accept- 
able region is outside the circle of intersection. Associated with each CARDIAD plot 
is a locus of centers plot, which indicates the centers and labels the frequency of 
each. Space limitations do not allow the locus of centers plots to be included with 
the CARDIAD plots in this paper - , but they will be mentioned and referenced as neces- 
sary. 

3. ILLUSTRATION 

Figs. 1 and 2 are CARDIAD Plots of a two-input, five-state, two-output model of a 
jet engine. The model is derived from a jet engine simulator called DYNGEN [3,4] 
and represents an F-100 turbofan jet engine with a fuel flow of 2.75 Lbm/scc. (full 


f 



throttle without afterburners). The inputs are fuel flow and exhaust area and the 
outputs are thrust and high turbine inlet temperature. This model is one of a 
series of such models presently being used in a set point study of an F-100 like 
jet engine. 

The analysis of CARDIAD plots proceeds as follows. Recall that, at any given fre- 
quency, the acceptable region is outside the circle if the circle is dashed or 
inside if the circle is solid. The first question of interest is whether the columns 
Of the system are dominant uncompensated. For this to be the case, the origin of 
the CARDIAD plot must be i.ncluded in all solid circles and excluded by all dashed 
circles, since the origin represents identity compensation of the column. This is 
not the case for either of the two CARDIAD plots of this system. The next question 
is whethe: the system can be made dominant by constant real precompensation. If 
this is the case, there will exist a point on the real axis which lies inside all 
solid circles and outside all dashed circles. Fig. 1 shows that the first column 
of the system can be made dominant at all frequencies by the choice of any constant 
x, which lies outside all the dashed circles of the CARDIAD plot. Fig. 2 shows 
that there exists no constant value that will make the second column of the system 
dominant at all frequencies. Thus, some form of frequency dependent precompensation 
will be necessary. 


& 


Before proceeding with dominating this system, some of the features of CARDIAD plots 
should be mentioned. One property is that a circle at a specific frequency in the 
plot for one column will be solid if the other column is dominant at that frequency 
and will be dashed if the other column is not dominant. From this fact it follows 
that the transition from one type of circle to the other in the CARDIAD plot for 
one column occurs when there is a change in dominance in the other column. Once 
again considering Figs. 1 and 2, these facts indicate that the second column is not 
dominant at any frequency since all of the circles in the CARDIAD plot for the first 
column are dashed and that the first column is dominant at low frequencies (until 
W-7) because the circles in the CARDIAD plot for the second colunn are solid for 
this and all lower frequencies. 


A second feature of the CARDIAD Plot is the effects of a column switch on the plots. 
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I 
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a 






that is, prcmultiplicatlon by a matrix with the only nonsero entries being off-dia- 
gonal l’s. The effects of such a switching of the inputs are that all solid circles 
become dashed circles, all dashed circles become solid, and the shapes of the column 
one and two plots are switched. The CARD I AD plots of the system with this type of 
compensation are given in Figs. 3 and 4. Note that the first column is now domi- 
nant at all frequencies without further compensation. This fact can be ascertained 
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Fig. 3. Column 1, GCs)**^ 


Fig. 4. Column 2, G(s)*K 


either from the fact that the origin in the CARDIAD plot for column one is included 
by all solid circles and excluded by all dashed circles, or from the fact that all 
of the circles in the CARDIAD plot for the second column are solid. 


Since switching the inputs makes one column dominant uncompensated, it seems a logi- 
cal first step in compensating for dominance at all frequencies. Thus, K is chosen 
to be fn i1 


•( 17 ) 


It is still necessary to make the second column of - the system dominant . From the 
CARDIAD plot for this column (Fig. 4), it is apparent that frequency dependent com- 
pensation will be necessary since there exists no point in the real axis which is 
included in all the solid circles of this plot. To design such’ a compensator, a 
function of s is fitted to the shape of the CARDIAD plot so that, at any given fre- 
quency, the compensator lies inside the solid circle associated with the same fre- 
quency in the CARDIAD plot. While it is possible to find a first order compensator 
that will make this column dominant, a second order compensator has been used be- 
cause this same compensator could also achieve dominance at four other set points 
of the model. K„(s) is the compensator that achieves dominance in the second column 
of e(s)%. 

1 r. -.742s - 9.59 1 


K 2 (s) 


.014s - o998s d- 1. 
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The CARDIAD plots of the system with this compensator are given in Figs. 5 and 6. 

It is obvious either from the fact that only solid circles appear in the plots or 
from the fact that all the solid circles include the origin that each column of the 
system is now dominant at all frequencies. 



Fig. 5. Column 1, G(s)*K(s) 


Fig. 6. Column 2, G(s)*IC(s) 


4. GENERALIZATION . 

The CARDIAD Plot approach to system dominance in the three- input, three-output case 

is similar to the approach In the two-input, two-output case. 

' 

The actual condition for dominance in the 3x3 case is the i column of a matrix 
Z(s) will be dominant If 

. kj^Cs)! > l jz^ool * * (20) 

35^i 

for all s on D. If both sides of this inequality are squared as. in the 2x2 case, 
then an equivalent condition is 

k<< s >I 2 >ri k,-< s >l [ • . . ox) 


Using inequality (21) , the condition for dominance in, say, the first column is 
Jz. (s){ 2 > |z„. (s)| 2 + Jz (s)j 2 * 2|z (s)||z (s)I. 


The cross term produced by squaring adds non-integral power terms to the dominance 





inequality for the 3x3 system. To circumvent this problem, the last term of ine- 
quality (22) is replaced by an upper bound. Since 


|z 21 (s)r + |z 31 (s) j / :> 2 jz 21 


(s) I |z 31 (s) 


(23) 


with equality when | z-. <s) 1 - jz-.(s)|, it is convenient to 'replace the last term 
of inequality (22) with the left member of inequality (23). This yields a suffi- 
cient condition for dominance. For column 1, the condition is 


| Zll (s)| 2 - 2lz^(s)l 2 - 2[z^(s)l 2 > 0; 


21 


31 


(24) 


1 = 1,2,3. 


(25) 


and the general form is 

1*^60 1 2 -2 l k 1 w| 2 »0. 

3=1 

From inequality (2A) , the derivation of the dominance equation for the 3x3 case 
proceeds analogously to the 2x2 derivation. The general form of the compensator 
used in the analysis is 


K(s) = 


where a 


ij 


1 « 12 (S) 

“l3 Cs) 

a 21 (s) 1 

S 23 (s> 

a 31 (s) a 32 (s> 

1 

+ y ij J - 



• (26) 


s=jtd 


Once again, the open loop transfer function matrix G(s) and the general form (26) of 
the compensator are evaluated at a specific frequency and multiplied to form Q(ju). 
Then, using inequality (25), a dominance inequality for each of the three columns of 
Q( 3 w) can be formed. For example, the first column of Q(j to) will be dominant at the 
frequency to if 


!q 11 | 2 -2| q „| 2 -2|q„| 2 > 0 


21 1 


*31' 


(27) 


and the dominance function for column 1 is 
, * 2 
^l^ X 2l 5 ^21 ,X 31* ^31/ " ^ ^ y* 


‘21 c 2 


2 , 2 2 
21 C 2 + X 31 C 3 + y 31 C 3 


* 2X ?,1 C 4 * 2y 21 C 5 + 2X 31 C 6 + 2y 31 C 7 + 2x 21 X 31 C 8 

* 2y 21 y 31 C 8 + 2x 2l y 31 C 9 " 2x 31 Y 21 C 9 > 0 


(28) 


where the constants c^-c^ are functions of G(s) evaluated at the frequency to. Simi- 
lar dominance functions can be derived for the other two columns. 

The maximum-minimum analysis is performed in two different ways. In the first 
approach, which will be referred to as the standard analysis, the variables of the 
dominance inequality are first paired by the entry in the compensator which they 
represent; and the maximum-minimum analysis is performed on each pair assuming that 
the other pair is zero. The resulting raaximums or minimums are 


- V c 2 ! 


r 21 


-c c / 


5 2 5 


X 21 


w The hessian for each pair of variables is diagonal and the second derivative test 
•once again reduces to a sign test. * 

1 The dominance analysis is repeated over a range of frequencies and CARDIAD plots 
* result. There is one plot for each off-diagonal entry in the compensator and each 
. entry is plotted assuming that the other off-diagonal entry in the column is zero. 

I Using CARDIAD plots generated by the standard analysis, dominance is achieved by 
setting one of the off-diagonal entries to zero while the other is chosen as was the 
case in the 2x2 design. 

9 There does not always exist a value in one off-diagonal entry of a column of the com- 
® pensator that will make the column of the system dominant when the other off-diagonal 
entry in that column of the compensator is zero. When this occurs, the maximum- 

I minimum analysis is performed by finding the full gradient of the dominance function. 
The hessian is no longer diagonal but the eigenvalues of the hessian are all nega- 
tive in Section 5, so the point that is found is a maximum. Design which is per- 
il formed on plots generated by the full gradient analysis involves both of the off- 
-fi diagonal entries of a column of the compensator, and functions must be fit to each 
to achieve dominance. 

1 A new symbol appears in the plots. At any given frequency, unless dominance can be 
achieved at that frequency with the other entry zero, a small triangle is drawn 
which shows the best that can be done towards achieving dominance. It should be 

f noted that the triangle can appear in plots generated by either analysis. In the 
standard analysis CARDIAD plots, if triangles appear in one plot for a column but 
not the other, dominance can be achieved by keeping the entry in which the triangles 

I appeared zero and using the other entry to achieve dominance. In the full gradient 
analysis plots, triangles appearing in both plots do not mean that dominance cannot 
be achieved. Given that one entry in the compensator is chosen exactly on the tri- 
ok angle .<t a -certain frequency, there is a radius of points around the triangle in 
■ the other plot that will achieve dominance; but since the size of the circle is a 
™ function of how well the other entry is fit to the triangles, such a circle could ‘ 
easily be misleading. Both of these points will be illustrated in the next section. 

5 5. THEME PROBLEM ANALYSIS 

f The following design is performed on the reduced order model of the theme problem 
with state feedback. The states being fed back are the two turbine speeds and the 
pressure P^. Dominance will be achieved using only precompensation. 

S The plots for the uncompensated system using the standard dominance analysis showed 
that the first two columns of the system could be made dominant with one off-diag- 
onal entry in each of the first two columns of the compensator zero. The third 

f column, however, could not be made dominant at any frequency with either one of the 
off-diagonal entries in the third column zero. Physically, this indicates that the 
principal effects of all three inputs (fuel flow, exhaust area, and guide vanes) 
are on the two speed states. To facilitate achieving dominance, a column switch was 
jg done by choosing the first compensator to be 

0 10 ' 

0 0 1 ., 

.10 0 . 

f . Figs. 7-12 are the CARDIAD plots of the system with this compensator and use the 
standard dominance analysis. The plots for the entries in the first column, Figs. 

7 and 8, show that the first column is dominant without further compensation, since 

f the origin of each plot is included inside all solid circles and excluded by all 
dashed circles. Figs. 9 and 10 are the CARDIAD plots for the second column. Fig. 

10, the plot for the 1,2 entry, has several triangles in it, indicating that, at 

I 






JJ.M 13.W 5 


Pig. 7. G(s)*K^» 2,1 Entry 


Fig, 8. 3,1 Entry 
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Fig. 9. GCs)*!^, 1,2 Entry 


Fig. 10. GCs)*^, 3,2 Entry 
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XI. GCs)*!^, 1,3 Entry 


diagonal entries in the compensator zero. 


Pig. 12. GCs)*^, 2,3 Entry 


Pigs. 13 and 14 are the plots for the third column using the full gradient rather 
than the standard analysis. The solid circles which appear at high frequencies in 
Pig. 14 are very important. Recall that the circle will only be drawn if dominance 
can he achieved while the other entry is zero. This means that by staying inside 
these solid circles, dominance can be achieved at the frequencies at which they 
occur while the 1,3 entry in the compensator is zero. Thus, ir. designing the 2,3 
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Pig i 13. G(s>*K t , 1,3 Entry (Gradient) Pig. X4* G(s)*K-, 2,3 Entry (Gradient) 






entry, the rat that is employed Is to follow the triangles; at low frequencies 
and ••y in idc ;> \ solid circles at the higher frequencies. 'If this is done, the 
design of the 1, r ntry will be simplified because it will only be necessary to fit 
the e:n ty to the • ow frequency triangles and have the function go to zero at higher 
Ci : ies. 

! ' ’ g this strategy, a lag compensator was designed to fit the 2,3 entry as 
.. . hed previously. The compensator entry that was chosen is 


V. 23 (s> 


-129.4s -1940.2 
.0365s + 1. 


At he same time, another lag compensator is fit to the solid circles in Fig. 9, 
U.u CARD I AD plot for the 1,2 entry. This was chosen to be 


k 12 Cs) " 


■ 0127 

,1162s + 1. 


Defining this compensator as K ? (s) with all the other off-diagonal entries zero, the 
c~ -rail compensation thus far Is K^Cs) = 


K 3 (s) = 0 


-129.4s - 1940.2 
.0365s + 1 


.0127 

•1162s + 1. 
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Figs. 15-20 are the CARDIAD plots of G(s)K 3 (s) using the standard dominance analysis 
The plots show that the first two columns of the system are dominant at all frequen- 
cies since in Figs. 15-18 the origin of each plot is contained by all solid circles 
and excluded by all dashed circles. Fig. 19 shows that the strategy applied in the 
design for the 2,3 entry was successful. To make the third column dominant, it is 
now only necessary to fit a compensator to the shape of the solid circles in Fig. 19 
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Fig'. 15. G(s)*K 3 (s), 2,1 Entry 


Fig. 16. G(s)*K 3 (S), 3,1 Entry 
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Fig. 21. G(s)*K(s), 1,3 Entry 

6. CONCLUSIONS 
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Fig. 22. G(s)*K(s), 2,3 Entry 
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The graphical CARDIAD method described in this paper has been effective on the Theme 
Problem. The authors’ experience indicates that it is an easily learned design aid 
which can be quite helpful in achieving dominance for realistic plants. A special 
advantage of the CARDIAD approach lies in the way in which it provides insight to 
the designer. The plots indicate whether or not it will be possible to achieve 
dominance with simple, lead-lag compensators. Examples up to this time suggest that, 
over the useful bandwidth, simple compensators are often successful in this regard. 
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It should be noted that this paper illustrates only compensator selection for domi- 
nance. Completion of the design is by classical means* For an example, see [5]. 

■■f particular note is the fact that compensator denominators having right half plane 
zeros do not necessarily lead to unstable controllers. This may also be seen in [5]. 
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Continued research on this class of graphical, interactive methods is in progress. 
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